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It is well known that an important class of problems in numerical analysis can be 
interpreted as “electrical’’ network problems [1, 2, 3]. These analogues have formed the 
basis of many analogue and digital methods of solution of certain partial and ordinary 
differential equations. The most extensive and diverse such interpretations have un- 
doubtedly been made by Kron [4~—9]. Furthermore, it has been pointed out by von 
Mises [10] that the network approach to numerical analysis offers promise for the numer- 
ical solution of physical problems for which classical methods, in spite of great effort, 
have been unavailing—problems for which future research demands a numerical answer. 
See also the recent work of Birkhoff and Diaz [11] utilizing network methods in non- 
linear problems. An excellent example of a successful attack by network methods is 
given by Oppenheim [12, 13]. See also the provocative work of Branin [14]. 

We shall here consider not the problem of devising topological models for physical 
problems, but the network problem per se. We begin with a resumé of its history. 

1. History of the network problem. Kirchhoff, as a student in Neumann’s seminar, 
made the first comprehensive investigation of the electrical network problem, publishing 
his results in 1845 and 1847 [15]. Kirchhoff proved the existence of a solution to the 
‘purely resistive’ network problem. In formulating the problem he was perhaps the 
first to write down the incidence relations of a 1-dimensional complex. Maxwell [16] 
pointed out (1865), however, that Kirchhoff’s formulation omitted the concept of 
potential. Maxwell remedied this omission and devised two methods for solving the 
network problem effectively, now termed Maxwell’s mesh- and node-pair-methods of 
solution. (The two methods may be described in Fig. 1 as that of inverting L’ and Y’.) 
In 1923 Weyl [17] presented another proof of Kirchhoff’s result giving deeper insight 
into its topological character. His 1932-33 lecture notes refer to currents and voltages 
as contravariant and covariant vectors. Eckmann’s 1945 paper [18] gave further insight 
into the problem. 

In 1933 and 1934 [19, 20] Kron utilized the tensor calculus to fit in the concept of 
“impedance” of stationary and non-stationary networks into the topological frame-work. 
This advanced the application of tensor calculus from field problems to circuit problems. 
Kron published a treatise on these developments for non-stationary networks in 1938 [21] 
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and for stationary network in 1939 [22]. See also [23]. Figure 1 is an interpretation of 
Kron’s view of the network problem; he had each mapping and viewed the problem in 
its general form. Actually the concepts of tearing and interconnecting appeared in these 
publications: in the process of writing down the ‘equations of performance” of a system. 
The idea of interconnecting solutions did not appear explicitly until 1953 [8]. Kron’s 
formulation of the network problem has been studied by many investigators [24-28]. 

An algebraic topological formulation of the network problem and a proof of the 
validity of the method of tearing were given by the author [29-32]. A condition called 
power definiteness [31] was shown to be sufficient for the existence of a solution to the 
network problem; it is shown in the present paper that the proof given in [31] actually 
utilized the weaker condition termed ‘‘ohmicness’”’. In fact, ohmicness (See. 2) turns 
out to be a necessary and sufficient condition for the existence and uniqueness of the 
generalized network problem, the problem actually encountered in tearing. 

Section 3 treats the question: what systems can be represented as a network problem? 
The answer is any system of linear equations can be so represented (but noé necessarily 
in a way suitable for efficient tearing.) In Sec. 4 tearing is briefly described. The in- 
variance of power under tearing is demonstrated in Sec. 5. Section 6 describes Kron’s 
“‘orthogonal’’ method of solution of networks. Section 7 describes the explicit construction 
of the solution of the network problem via tearing. Section 8 describes A-partitioning: 
this is a way due to Kron of solving a system of linear equations which in general requires 
fewer operations than the customary inversion of a matrix. Sections 9, 10, and 11 give 
explicit solutions to three network problems, counting the required number of operat ions 
for tearing, K-partitioning and standard inversion methods. 

2. Necessary and sufficient conditions for the existence of a solution to Kron’s 


problem. We give a more general formulation of the network problem than in [31]. 
A network is defined as a 1-dimensional regular complex. Let K then be a network and 


let L be a transformation of its group C’(K) of 1-chains into C,(K) its group of 1-cochains, 
over the complex coefficient field. Then K and Z define the diagram shown in Fig. 1. 
Here H'(K) is the first homology group, P’(K) is the group of 0-boundaries, the image 
of C'(K) under the boundary operator; H,(K) is the first cohomology group and P,(K) 
the group of 0-cochains modulo the kernel of the coboundary operator 6. l'urther C is 
the natural homomorphism of H*(K) into C'(K) [C'(K) may be considered as H'(K, K°), 
where K” is the 0-skeleton], while C, is the dual transformation of C,(A) into //,(K); 
d@ denotes the boundary operator (strictly speaking, the boundary operator followed by 
the projection of C°(K) into P®(K)). Finally L’ = C,LC, Y = L™' and Y’ = dY6. The 
network problem then is: Given arbitrary elements J’ of P"(K) and e’ of H/,(K) to find 
elements J of C’(K) and V of C,(K) such that 0/ = 1’, C,V = e’ and V = L,/, i.e. to 
find a general expression of V and J in terms of arbitrary e’ and J’ and the transformation 
of the above diagram. A pair (J, V) satisfying these conditions for given J’, e’ is termed 
a solution. In [31] it was shown that a solution always exists and is unique if L is power 
definite. Professor N. E. Steenrod pointed out, however, that the proof actually utilized 
only the following weaker condition: 

Definition. Let V be a vector space of complex dimension n and V, , the dual 
space of linear functions on V. A linear transformation of V into V, is termed ohmic 
if for each non-zero v in V, (Lv), « (v) ¥ O, i.e. L is ohmic if it never maps a non-zero 
vector v into a function which annihilates v. This condition is more general than power 


definiteness, which in turn is a more general condition than positive definiteness. 
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Professors R. J. Walker and T. 8. Motzkin of Cornell University and UCLA have 


suggested the following interpretation of ohmicness: Let L = M + i where M = 
(L + L*)/2 and N = (L — L*)/2i are hermitian. Splitting these into their real and 
imaginary parts MZ = M, + iM,,N = N, + tN, we set 
aa (i -M) pa [Ns — Ma), 
M, M,) N. Ni) 


Then L is ohmic if and only if the two quadrie hypersurfaces 
z,Ar = 0 «.Be = 0 


have no real intersection. 

The introduction of an isomorphism L between a vector space V and its dual V, 
is equivalent to the introduction of an inner product in V for if v, w « V, we may define 
v-w as (Lv)-w. If L is ohmic this inner product is not necessarily commutative. If V is 
a vector space over the reals and L is positive definite then the usual properties of an 
inner product are enjoyed. Although the results below guarantee that tearing is valid 
under the more general condition of ohmicness, all the applications to date have been, 
to the author’s knowledge, for L positive definite. In Kron’s model of Schrédinger’s 


equation, however, L is ohmie but not positive definite. 


Theorem 1. If L is ohmic, the network problem has one and only one solution. Yor 
let 7’ ¢ H'(K) and L’v’! = 0; then 2’-L'7’ = 0. Letting 7 = C’d’ we see 7’- Lr’ = 1-Lis 


i = O and since C is an isomorphism 7’ must also be zero so that the kernel of L’ is zero 
and L’ is thus 1-1. Similarly Y’ is shown to have an inverse and J = C(L’)'e’ + 
Y6(Y’) ‘J’ and J LJ will constitute a unique solution. Uniqueness follows by the 
same argument as in [31]. 


IXron’s approach to networks (whether stationary or not) is always from the following 
point of view: He considers not just a single network with its impedance tensor, but the 
set of all possible networks into which it can be transformed through the group of tearing 

and intereconnecting—transformations. His method of tearing consists in deducing 
from the “equations of state” or the ‘‘equations of solution” for one element of the set 
(usually a simple one) the ‘‘equations of state” or “of solution” for any other element in 
the set by a routine procedure governed by the induced transformations. Stated another 
way, a given system is torn apart, solved in pieces, then the piecewise solutions are 
interconnected to yield the solution for the whole. 

We wish to find a condition on Z such that a solution will exist for all networks K 
compatible with L (we say K is compatible with L if the dimension of C'(K) is the same 
as the dimension of 1); i.e. regardless of how the network is hooked together we still 
wish to be guaranteed the existence of a solution: we term this the generalized network 
problem of Kron. 


Theorem 2. A necessary and sufficient condition that the generalized network problem 
have a solution ts that L be ohmic. 

Theorem 1 establishes the sufficiency; for necessity, let v be a non-zero 1-chain such 
that v-Lv = 0; now construct a network where v is a cycle, then it is clear that L’ will 


not be an isomorphism. 
3. What systems can be represented as a network problem? ‘This question is of 
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especial interest to those concerned with explicit numerical solutions of their problems 
(see [33]). If tearing is a shortcut method for solving the network problem, can it be 
used to solve arbitrary non-singular systems of linear equations? The answer is that 
there is an infinity of network problems which can represent a given such system but 
there is no a priori method for picking out a ‘‘good’”’ such representation. We give here 
one such representation (which does no good whatsoever from the viewpoint of com- 
putations). Let Lx = b be a non-singular transformation (Z can be a “function” of 
xz and b), L of dimension n. Let K be a network in which the natural mapping C: H'(K) > 
C'(K) is an isomorphism; we then consider L as a transformation from C'(K) into 
C,(K) and thus C,LC = L’ is also an isomorphism. (K consists of n 1-spheres with 
exactly a single point in common, a “‘bouquet”’ of 1-spheres.) By Theorem 1 the network 
room”’ to tear (see Sec. 4). 


‘é 


problem automatically has a solution, but now we have no 
We must invert L’ which is the same as inverting Z. Any non-singular system of finite 
difference equations can thus be interpreted as a network problem, which harmonizes 
with von Mises’ use of the term. The more usual mode of viewing this is the “all” june- 
tion-pair of equivalent network. This result is related to that of Lyusternik [2] and Bode 
[34] on the question of realizability. 

In Kron’s approach one attempts to represent not a given set of equations but rather 
a physical system by means of a topological model. While the topological model can be 
made to produce a purely algebraic system, the process is not reversible: the topological 
model conte’1s more structure, more information, than a purely algebraic system. It 
has been recently suggested by Professor Synge (Princeton University Colloquium 
Lecture, Feb. 7, 1956) that certain fundamental problems in quantum mechanics contain 
inherent topological structure which must be understood “or their solution. 

While the all-resistance or resistance-capacitance models of Poisson’s equations are 
widely used, not so well known are Kron’s more subtle topological models for a large 
class of partial differential equations of mathematical physics (both classical and quantum 
physics) such as the field equations of Maxwell; the wave-equation of Schrédinger; the 
Navier-Stokes equations of compressible, viscous flow; the elastic field, etc., [3-9], 
[35], [36]. 

4. Tearing. Let K and K”® be regular l-complexes and ¢ a continuous map of K”* 
into K which induces an isomorphism Q of C'(K*) on C’(K). Figure 1 of [32] depicts 
the maps induced by ¢ on the reduced homology and cohomology sequences of K and 
K*. If L is an ohmic isomorphism between C’(K) and C,(K) then L*’ = Q,LQ defines 
an ohmic transformation of C’(K*) on C,(K*). Thus if K is compatible with the tensor 
L so is K® and thus associated with L is the collection of all K* for which there exists 
a transformation g inducing an isomorphism of the l-chain groups. Note that ¢ is not 
necessarily a topological map. In [32] a description is given of how a network problem 
in K defines a network problem in K* and how computation of the solution for K” allows 
for the deduction of the solution for K. 

First we shall discuss an ancillary question—the invariance of the power under 
tearing and interconnection—then we shall give an explicit description of the construction 
of a solution via tearing and conclude with some concrete examples where tearing gives 
quicker results. 

In recent unpublished papers Kron has solved eigen-value problems of large-scale 
physical systems and linear programming problems by means of tearing. He has likewise 
utilized tearing to obtain analytic solutions of certain large-scale problems, e.g. involving 
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resonant frequencies of nuclear reactors. Indeed the savings in straightforward numcrical 
solutions seem to be trivial compared to these developments [37]. 

5. The invariance of the power under tearing [33]. Let J « C’(K) and V = LJ. 
The power is defined as the product V-J. Explicitly, V-J in matrix notation is V*-J = 
J*L*J; (power as here defined is the sum of the disstpative (V,J* + V*J)/2 = J*(L + 
L)*J and reactive (V,J* — V*J)/2i = J*(L — L*)/2¢ power, the latter multiplied by 2). 

Theorem. Power is invariant under tearing and interconnection. For let g be a map 
of network K* on network K which induces an isomorphism Q of C'(K*) on C'(K). 


The induced transformation 

C'(K") * CK) 
on the group of 1-cochains is then 

C,(K*) & C,(K). 


Let J* be an element of C'(K*): the power corresponding to this element is V*-J*. 
But if J/* = Q/ and V = Q,V”* then the corresponding power in K is V-J = (Q,V")- 
(Q’J°) = V*-J*, aed. 

This theorem settles a problem which apparently long has plagued the literature 
in network theory, namely the invariance of the power under the transformations 
induced by interconnection. The concept of tearing as a form of .transformation has 
been emphasized by Kron [36]. Still many investigators, including Kron, have attempted 
to deduce the law of transformation from the invariance of the power. 

The fundamental geometrical fact about tearing and interconnecting is that there is a 
one-to-one correspondence between the branches of the two networks: the transformation ¢ 
of inlerconneciing and its “inverse’’, tearing, explicitly establish this correspondence. If 
the complexes are simplicial, we may restrict the transformation of interconnection to 
be a simplicial map: but the really essential feature is that it establishes this 1-1 corre- 
spondence between branches: expressed another way, the transformation ¢ shall induce 
an tsomorphism of the 1-chain groups. In algebraic topology the laws of transformation 
of 1-chains and 1-cochains (currents and voltages) are well known and completely 
understood (see [38]): the underlying geometrical transformation dictates the trans- 
formations (P, Q, R, P, , Q, , R, of [32]) in the algebraic superstructure. 

And thus the invariance of the power is a consequence of these transformations and 
not a cause of them. The other approach therefore was forced somehow to smuggle in 
the assumption of the invariance of the power. For instance, le Corbeiller [28] in his 
fine little book on Kron’s method employs very detailed but circuitous reasoning in this 
attempt. 

6. Kron’s ‘‘Orthogonal” method of solution of networks. Preparatory to the 
construction of an explicit solution via tearing we must first consider Kron’s so-called 
“orthogonal’’ method of solving networks. Maxwell’s mesh method and node-pair 


c F) 
Oo —— > HK) ——*> c!(K), ——— p® (k) —e 0 


Cy § 
0 +— H, (Kk) — Ck) ————P, (k) @— 0 


Fia. 1. 
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method of solving networks in modern terminology amounts to inverting L’ or Y’ 
respectively (see Fig. 1); the method we describe here stays within the spaces of C'(K) 
and C,(K) and although in general this method requires exactly as much work as the 
mesh or node-pair method it gives a new twist to the subject, to pave the way for tearing. 

It may be seen that the following method gives a solution. Let C’(K) be written as 
the direct sum of its cycle group R' plus some other group, S', chosen at one’s pleasure: 


C(K) =R'+S'. 
Let a basis be selected for R' and S', so that an element J of C' can be written in coordi- 


nate form 


J = Me 
I) 
i, I being the coordinates of J in R', S‘. Let us choose the ‘‘same”’ basis for C,() (see 
[32], also Cartan’s cosimplex, [39], p. 39): this choice of basis effects a splitting of C,(A) 
say 
C.0n) = R, + S, 
(of the four groups R', S', R, , 8S; , only R’ is natural 
Let an element V of C,(K) be written therefore 
v=|°| 
LE J 

with e, EZ, elements of R, and S, , or rather e, FZ, being the coordinates of V in Rk, , 8, 
for the bases chosen. Note that R,; , which corresponds to the dual basis for R’, is disjoint 
from the subgroup of coboundaries and that in general S, is not the coboundary group. 
Nevertheless the natural homorphism of R, into H/,(K) is an isomorphism. (The reader 
should now refer the argument to Fig. 1.) Thus an element e’ of /,(K) may be pulled 
back to an element e of R, . Likewise an element J’ of P’(K) can be pulled back under 
to an element J of S'. Thus we have elements e, 7, such that C,e = e’ and d J= I’. 
(Note that in some actual problems such an e is frequently given in the original specifica- 
tion of the problem). We shall write the ohmic transformation between C,(K) and 
C'(K) in the form 

| 2] iY, Ys; 
r 


LY; Yq 


sj * (6.1) 
E 


) 
| | 
| 
| 
J 


We shall call this the first normal form. The solution may now be obtained by inverting 
Y, (Y, is ohmic by heredity) 
E = Y,'(I — Y~,¢). (6.2) 
Then i may be obtained as Y,e + Y.E. The above formula for F will be termed the 
“orthogonal” form of the solution. This is the form of the solution we will use for the 
torn system in the next section. 
Alternatively if we write the ohmic transformation in the form inverse to 6.1, termed 
the second normal form, 


| °) = in Ly) 0) (6.3) 
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This is the form assumed by the problem for the interconnected system after interconnec- 
tion [see Eq. (7.5)]. In this form a solution is effected by inverting L, : then 


i= L;'e — LJ), 
K= Lyi + Ll. 


7. Explicit construction of solution of network problem via tearing. An important 
consideration not previously discussed is a convenient method of construction of the 
solution via tearing. We shall be working with coordinates which will transform then in 
the direction opposite to the transformations of the vector spaces themselves. Let K* 
be a torn version of K so that there exists a regular map ¢ (simplicial map if complexes 
K and K”® are simplicial) inducing an isomorphism ¢” of 1l-chain groups. 

We shall split C' = C'(K) and C} = C'(K") into the following direct sums 


C; =M,+P,, 

C'=M°+M'+P. 
Hlere M, is the subspace of cycles of Cj ; M° is the isomorphic image of 1/, under the 
transformation ¢* induced by ¢, M' is a subspace of cycles whose direct sum with M° 


“fills out”’ the subspace of eyeles. P, and P likewise fill out Cj} and C* respectively, but 
are so chosen that g*(M' + P) = P, , and thus ¢* splits in the following fashion 


M°’|M' P 


a ee U’ = identity 


[t is clear that such choices can be made. 
Let bases be chosen for each of these subspaces. Then in coordinate form, for zt , [s 
elements of M, , Py and 7°, 7, J elements of M°, M', P, ¢* takes the explicit form 


is) _ (U 9) e 
7 cars (7.1) 
I} [0/1 Q\ |, 


Let the “same” bases be chosen for C, = C,(K) and Cf = C,(K"*) inducing thus a splitting 
of each of these spaces: C, = M, + M, + P, , C% = M* + P*; let e ,e,, Eande’, 
E* denote elements of these subspaces in the order implied. Then the induced isomorphism 
g, on the l-cochain spaces is written 

eo| U|O|{e* 

e| | | || |. U = identity (7.2) 


} E} O Q, E « 


Suppose now that the K* network problem is written Kron’s “orthogonal” fashion, 
with respect to the basis chosen for K* 

. | rs * 

a } } e - 

4 ee ek (7.3) 
, ” 
ALG. 


‘ 
nn * 


~ 
* 
~—_ 
‘ 
~ @ 
~— 
oa 
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This is the first form for the “orthogonal” system [see Eq. (6.1)]. The solution is obtained 
by inverting Y4. Let L’ = (Y4)"' and /* = J, — Y%e*. Then 

E* = L*I’. (7.4) 
(For the purposes of interconnecting, however, it is, of course, quite immaterial what 
manner is employed to compute E” in the form L*I*). Equation (7.4) represents a solution 


of the K* network problem: we now construct the solution of that of K via the interconnection 
transformations. Utilizing successively (7.2), (7.4) and (7.1), 


fe _ Q.E" = Q,L*I* = Q,L*Q : . (7.5) 
LE I 
We now split Q,L’Q into its component parts 

: Vf 
¢ - ps _ v7] ; (7.6) 
LE \Ls L, L F 


This is the second form for the ‘orthogonal’ system [see Eq. (6.3)]. The work that 
5S A t 1 

remains to be done therefore, after the interconnection, is the inversion of LZ, , which 

may be termed the elimination of the newly created cycles. Thus we write as the final 





solution 


Il 


i = L;'e, — LD, (7.7) 
E=L,yi' + Ll. (7.8) 


Of course (7.7) may be substituted directly into (7.8) to obtain FE = (L, — L,L;'L.)I + 
L;L;7*e, but from a computational viewpoint it is advisable not to do so. In the next 
section we shall consider in detail the question of the number of operations performed. 
Reviewing the actual steps taken for solving by means of tearing we see that the following 
operations had to be performed. 


I 


A) In (7.3) Y 4 was inverted: (Y4)~* = L’. 

B) In (7.5) the multiplications Q,L*Q were performed; it is to be recalled, however, 
that Q is a matrix with elements 0, 1 or —1 so that no arithmetic multiplications, 
only additions were performed. Furthermore (see below) it actually would be 
necessary only to form the Z,-part of this product. 

C) In (7.7) L, was inverted. (Actually one need only ‘‘solve” Z, , say in the manner 
of K-partitioning below). 


Summing up then, tearing involves inverting two matrices Y4 and L, . Now let n denote 
the number of branches of K (which equals the number of branches of K*). Let p* be 
the dimension of H'(K), the first Betti number of K, and p, the first Betti number of 
H*(K*). The dimension of Y{4 is n — py and that of L, is p’ — pg : if K were inverted 
by the Maxwell node-pair method one would have inverted a transformation of the 
dimension n — p' which is less than n — py . How then can it be claimed that it is possible 
that the method of tearing may require fewer arithmetic operations than the standard 
methods? The answer is that tearing must be done in such a way that the Y4, to be 
inverted split in some convenient fashion. This depends upon the topology of K* as 
well as on the ohmic tensor. In Secs. 9, 10, 11 we shall exhibit cases where it does pay 


off to tear. 











. 
. 
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8. K-partitioning. We are about to compare in detail the number of operations 
required to solve a few selected network problems by standard methods, by standard 
partitioning, by an algebraic form of solving equations suggested by tearing, and by 
tearing. We shall therefore preface these examples by a few remarks on partitioning 
and obtaining the general solution to a system of linear equations. See [40]. 

Let 

Mz =c (8.1) 


by the matrix form of a non-singular system of n equations in n unknowns. The usual 
manner of obtaining the general solution is to compute the inverse of M by one method 
or another, writing the solution in the form z = M~‘c. Now to invert M it may be 
expedient to partition M, but let us instead partition the whole system (8.2) into the form 


) 





A B 
A Bi\\u _ |@ (8.2) 
C Dj Lb 
where 
A Bl _y iy — fl =e. (8.3) 
c D} Lv b} 
An easy computation shows that the general solution may then be written as 
_= b — CA“a) a 
U R(b CA a), p= (D —CA 'B) : (8.4) 


A-‘(a — Bo), 


(8.4) is a general solution in the following sense: for arbitrary a and b, v is computed 
according to the first formula; it is then substituted into the second to give u. This is 
just as good a form of the general solution as (8.2). If (8.4) is accepted as a legitimate 
form for the general solution, we see that the only multiplications involved are in the 
computation of A~‘, the formation of CA~'B, and finding the inverse of (D — CA~’B). 
This form of the general solution requires in general significantly less work than finding 
the straight inverse. Incidentally, if the transformations are multiplied out we obtain 
the standard formula for the inverse by partitioning (U = identity matrix): 


_ {A(U + BRCA™") —BR|{a 8.5) 
Lv] | —RCA™ R J\b 


The form (8.4) is a purely algebraic trick, simple but useful, suggested by the method 
of tearing. We shall term this Kron’s method of partitioning or K-partitioning, and will 
compare it with standard partitioning and with tearing in solving network problems. 
The disadvantage of the component solutions found by K-partitioning, compared 
with those found by tearing, is that the former cannot be re-utilized to interconnect 
the solved system with other systems to form the solution of “‘supersystems’’. No top- 
ology enters into partitioning a matrix. On the other hand, the component solutions 
established by tearing represent actual physical systems and consequently all partial 
inverses can be re-utilized to build up factorized solutions for still larger systems. 
Kron has used such sequences of partial structures for pyramiding not only numerical 
solutions, but also eigen-value of analytical solutions, as well as problems in optimization 


Uu 


| u | 


(largely unpublished). 
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We shall use throughout von Neumann’s rule-of-thumb, that it requires n* multipli- 
cations to invert a matrix of order n. 

9. Tearing versus partitioning: first case history [41]. Consider the network in 
Fig. 2. Here n = 6, p' = 3. We assume that L is purely resistive, and thus that it can 


S 














Fig. 2. Network K. 


assume diagonal form. Solution of this by either the mesh or node-pair method requires 
the inversion of a 3 X 3 matrix of the form 


iz O a] 
lO ZS £ ‘x’? means a non-zero element. 
2 £ £ 


According to our adopted rule-of-thumb, this requires 27 multiplications. My computa- 
tions yield the following results: 


Number of multiplications required to solve problem 
I Y yi 


























Tearing 3 
K-partitioning S 
Standard partitioning 22 
Standard inversion 27 
a7 
as” 
lg 24 
\ ov 9 
37 
_ 
K 
iw * 
' 2 
2'% 
“ 
K 


Fig. 3. K and K* with bases selected. 
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Solution of K by tearing. We shall tear K into K* as shown on the top of Fig. 3. 
In the figure for K*, 1, for example denotes a “current”, a cycle impinging on the two 
branches which make up the cycle; 1* indicates a current passing through only the single 
branch indicated, in the direction shown: this is a 1-chain consisting of this single directed 
simplex. Similar meanings are attached to the other labels. 

First we must solve the K* network in the first orthogonal form (6.2). A most econom- 
ical method of doing this is to start with the completely torn system K* below and 
interconnect from K* to K*, for L must be given a definition completely in terms of 
coordinates, and it customarily is given for the K* system, Fig. 4. 














a® 
, 5 \* a e* 
3% 
Fig. 4. The completely torn network K*. 
Suppose Y* = 17’ assumes the form 
[Y. 0 
y* ai Ya 
if 
eo 
| 0 Yo 
In the AK* system the ohmic transformation is written 
J* = y*y* 
where J* = {J'*, --- , J*}, Vy = {Vi, , --- , Ve,} the components related in the 


obvious way to the symbols of Fig. 4. (J* and V* do not denote complex conjugates 


in this section.) 
Let y be the simplicial map of K* on K*, whose induced isomorphism on the groups 


of 1-chains is 


1* 2* 3* 4* 5* 6* 
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26 





Matrix of interconnection from K* to K*. 
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We may write this in the form 


-# 
4 


{7s 

{J 

= P,J*, 

where {7*, 7} = {7'*, 1°*, 1°*, I**;7'*7"} and P, is the matrix shown above. The inverse 


of these Q’s may be written by inspection, for if 


cr 


p= (0 7 
Ge 
P, arbitrary, the U’s of different size in general it is immediate that 
Pp" a U —P, 
LO U 


Kron has used P;' instead of P, , expressing ‘‘old currents in terms of new’’ and he 


writes it down by inspection of the network; see interconnection from K* to K below. 


Then 


and so 

s 3 
: = P,J* = P,Y*V* = pvep|* ; 
1 e* 
Written more fully, 


’ 


"|= [¥ Y*%)(E* 
(? ly’ Y%Jle® 


E* = (Y‘)"(I* — Y%e"). 


The inversion of Y* involves two operations: Yj, , Yo, , for Y3: and Y;: are assumed 
known. (It is assumed that Z and Y are known for the completely torn system.) Let 


| Lis O 
ye nte | le 
‘. 
| O ha! 
and I? — Y¥*%e* = I, ; we have 
vil (9.1) 


= ° ° . . . re , 
Now we are prepared for the interconnection transformation between K" and K. In 


contradistinction to the interconnection between K* and K”’, we write the “old currents 


in terms of the new”’, i.e. in the form (7.1). Referring to Fig. (3), we have 
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Matrix of interconnection from K¥* to K. 

















[7 
(7 Q | O} |i° 
")- Ie (9.2) 
i” O | U) {a 
. 
and the corresponding transformation on the cochains is 
(E 
| 
eéo|  (Q. | O| {ER 
—| = a = (9.3) 
iO | U) le* 
2) 
Thus from (9.3), (9.1) and (9.2) 
(RF) 
a = Q,E* = Q,L"I* = Q.1"0! ‘). 
Leo) Le) 
Now 
fe | Lis] 
Q,.L°Q = | Los oa * * 
Les | Lss ibe L,} 
| _—_————+ 
Lis Lee Lae | Ln 
where 


L, _ Lis + Los + Le + Lis ; 


the formation of Q,L*Q requires zero multiplications) or 
1 I 


ote leh 


rc 


BE L, 


L; 





“0 
Lo L, t 
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Remaining is to invert L, : this takes one division = one multiplication. Now as with 
(8.4) we write the general solution in the ‘“factorized’”’ form 

= Li'(eo — L,I), 

E=L,I + Li. 


I 


ll 


This completes the description of the solution of the K network problem. The method 
of tearing K utilized here is quite distinct from Kron’s method, several lengthy steps 
are eliminated; both depend, however, on the same fundamental topological principles 

One further remark can be made: as Kron has pointed out, the full affair Q,L"¢ 
need never be formed; only the piece L, , the lower right corner, need be formed. In 
fact, all that need be stored is: 


? 
Q 


4 


L,— 1. 


10. Tearing versus partitioning: second example, Poisson-type network. ‘This is a 
Poisson-type network except that, more generally, we assume that the “resistances” 
in each branch are different; Figure 10.0 shows the network K above with an indicated 
basis and a network Ky, below into which K is torn, also with an indicated basis. We 
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Fig. 5. Left—K and basis for C(K). Right—K* and basis for C(K* ). 








assume that LZ is diagonal for the completely torn network. First we give complete 
details on our solution of K by tearing. 


10.1. Solution of problem by tearing. 


10.1A. Solution of the K* network. It takes fewer operations to solve this network 
in the second normal form, i.e. we will solve for the “mesh currents’’, since the principle 
pieces are of identical character (note carefully that we have not assumed that J is 


the same for each) we shall compute the number of multiplications for one and then 
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# 4* 5* 


(1) 
x! (1) 
+ Ke 


Fig. 6. Upper left rectangle K¥'’’ of Ky and K{” , K{? completely torn. 


multiply the total by 4. To write down the appropriate equations we consider the upper 


left rectangle Ky’ of K, let K, be this network completely torn up. Then the matrix 


r vl neeti ) *(1) > r(W) 5g 
of interconnection Ky A, if 








a’ B hy - 3° .° 5° 
a*| 1 | | | | | | 
; es meee eee ies pis 
ev) {fy | | } ft | 
mf af of af of ff | 
= = _| = 7 im > — i —_— 

2* | —] | 1 | | J | | | | =Q 
See ae ae a 
| “ 
rim} | | tt | 
ot teak ft £ Wek 1 
ia ed 
| | 1] | | | 1 | 





and if L, for AY” is diagonal with elements 
Ly = {Loe ; Lge ’ | ; Lew ’ Lge ’ Lge ; Le} . 


Then letting L., + Li, + La, + La, = Las and Lg, + L., + Ls, + Ls, = Lops , we 











res \s 
[ Las —Loe | —Ly Ly —Ly ; 
ie. ite a 
LS? = Q,L°Q = —Ly | Ly 
th | % 
—Lye —Ly Kes 
—L,. | | 
Ly} « 
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Let L‘” be written 





Tra) | (1) 

#1 | #2 

is =. isi, ks 

4s — | 1 

(1) a) 
“#83 | Reidy 

The second normal form 

f ,@)) fra) (1) ( 
# | | Lisi 4#2|| Us 


| 


— 
| 





a) L (1) (1) (1) 
LEs J Liss Leds 
requires first the inversion of LY). For simplicity (and with no loss of generality) we 
assume ef? = 0. 

Then 


1(1) (1) (1) 1)\-ly (1) 1 4 (1) 
E, = [L < L e3(L 1) L salls = Lay #- 


A tally of operations yields: 


Operation Multiplications required 
Inversion of Lg, 8 
Formation of Ly7'L. 12 
Formation of L,(L;*L.) 30 
Total to solve K{” 50 


Now assume that we do this for each of the four main pieces of K* : the total number 
of multiplications is 4-50 = 200 to solve K, . The final solution for Ky is a direct sum 
of the solutions for each piece and the remaining disjoint pieces, the “pigtails”: for 
convenience this latter is split into two pieces L* relating to branches labelled 21 through 
24 and L* for 25 through 32, L* and L* both being diagonal. 
































wey Lee i fe 
E? LS, 1 
E? | _ L*, 0 lrg 
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Further let 
7 a i i a i - 
Ly _ Lay) T+ Lis) + L; 
and 
# + td 
M = L\ ot. L, ° 
This is the solution for K, . Now we interconnect to get the solution for K. 
# 
10.1B. The interconnection K, — K. Referring to Fig. 5 we see that the pertinent 
part of the interconnection transformation from C(K,) to C(K) is (we omit the part 
injecting the cycles of K, into the cycles of K—it would correspond to an identity 
matrix): 
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Interconnection tensor 7’: CK* — CK 


Going through the usual routine we perform the rearrangements T,I,T (this trans- 
formation involves at most addition, and actually need not be completely carried out— 
only the lower corner—the rest can remain in factorized form, a la Kron, see Sec. 8 


above). 
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[ Li LR) 

L=T7,M,T = |———— N = R,LiR + S,L*S. 
(LiR), N 

Again it will be observed that the formation of L involves no multiplication; N is the 

only “piece” of the transformation LZ that we must compute explicitly anyhow. 

The equation for K thus reads 

fe) (Li LR 1 

ia» 

Lo) \(L’R), NJ te 


The final remaining step is the “‘solution’”’ of N. 


10.1C. The elimination of the newly crealed mesh currents or the “‘solution” of N. 


Now N is of the form 














l 
y x x x x 
bic x =u Ss x 
| 
| 
c|x x x x|x 
| ae : 
d | : ££ 21x = Mi. 4 
\ Prmenrens San rm piled Bile 
e %Riz xs xz Nz | Naz 
| 
f Sie & = & 
| 
g|x ‘se = we 
h|x x | xX xX x | 








which we partition as indicated. Again, we do not compute the inverse of N but we 
compute the factorized form of the solution (Sec. 8). We tally the results, the method 
already having been explained by earlier detailed examples. 


Operation M ultiplications required 


Inversion of N, 56 
N;N;'N, 48 
Inversion of (V, — N;Ny7'N.) 64 


T otal number of multiplications to solve N in factorized form /168| (as compared with 
512 if N were inverted by standard methods). 


10.1D. Final tally of number of operations required by the method of tearing here em- 


ployed. 
A. Solution of Ky, network 200 
C. Solution of N, the elimination of the 
cycles created by the interconnection 168 


Total number of multiplications required by 


tearing solution 368 
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Next we discuss briefly a partitioning type of solution of K. 
Let us now refer to the original network K. 


10.2. Solution of K by K-partitioning. 
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We shall consider it in its mesh-equivalent form, that is, we “invert” the isomorphism 
(Fig. 1) of H'(K) “ H,(K), since L’ is 16 X 16 whereas Y’ is 24 X 24. Of many bases 
tried, the following checkerboard arrangement seemed to be the best. 
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Fia. 7. 


Checker board basis for H'(K). 


The transformation L’ for this basis assumes the following form, which is partitioned 


as indicated 
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Fia. 8. L’ for checkerboard basis. 
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“‘y”’ indicates the 
presence of a non- 
zero element, 
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Here is the tally of operations: 


Operation Multiplications required 
Aq 9 
CAB 98 
(D — CAB)” 512 [(D — CA™'B) is full] 


‘otal number o iplications required solving as 
Total ber of multiplications req 1 solving 


mesh equivalent network using K-partitioning (618) 
Next we see how many additional operations would be required if the actual inverse 
of L’ were computed by standard partitioning. 
10.3. Solution of K by standard partitioning. Now we must add to the above total 
of 618 the number of multiplications needed to form the inverse matrix: [see (8.5)]. 


Multiplications required to reconstruct actual inverse of L’. 


Operation Multiplications required 

R(CA™*) 8(24) = 192 [(C.A~*) has 24 non-zero elements] 
B(RCA™) 192 

A(U + BRCA™") §*° = 512 

BR 192 


Total 1079 


Thus for standard partitioning the number of multiplications required is 618 + 1079 = 
1697} = number of multiplications required by standard partitioning. 
10.4. Finally if one simply solved K by the most naive method, by inverting L’ by 
standard methods, then 16° = 4096 multiplications are required 
Final Tally 


M ultiplications required for solving K by various methods. 


Tearing 368 
K-partitioning 618 
Standard partitioning 1697 
Standard inversion methods 4096 


It will be remarked that the ohmic transformation L in this example was not assumed 
to act identically on the four principal pieces of K, (see Fig. 5) and thus the solutions 
for each of the pieces had to be solved separately. In the next section we shall consider 
the case when L does act identically on the four separate pieces. It will turn out that the 
method of tearing can take better advantage of these symmetries than can partitioning, 
even K-partitioning; tearing being a topological method can better accommodate to a 
topological situation than partitioning, a purely algebraic device. 

11. Tearing versus partitioning: third case history. Assume now for the network 
problem of Sec. 10, that Z acts in an identical fashion on the four principal pieces of 
K, . If we considered K as an electrical network, then we would assume that the imped- 
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ances of corresponding branches in these pieces are identical. We shall enquire as to 
how each of the methods can profit from this symmetry. 

11.1. Tearing. Referring to 10A, to solve the piece K§” required 50 multiplications, 
we make four copies of this; then we will assume that any degree of symmetry inherited 
by N from the symmetry cannot be taken advantage of, so that, as in 10.1C, the inver- 
sion of N requires 168 multiplications. 


Grand total for tearing = |218 


11.2. Solution by K-partitioning. A. For the checkerboard basis. Here unfortunately 
the checkerboard basis, Fig. 7, cannot take especial advantage of this symmetry in L: 
the inversion of A, Fig. 8, requires 2 instead of 8 multiplications and otherwise the 
counts are the same so that for this new problem, checkerboard basis with K-partitioning: 
613 


B. For another basis. Let us choose the basis shown below 








1}; 2/14] 3 





13; 19| 10; 4 





8} 12{ 11) 15 





7/16) 6] 5 




















Fic. 9. A basis for partitioning. 


Then L’ assumes the form shown in Fig. 10, which we partition as indicated. 
' 


Count of multiplications by K-partitioning 


££ - = 8 
CA" 24 
(CA~")B 33 
(D — CA™'B)”: 8° = §12 


Total [577| 


If L were not the same on each of the Kj then the inversion of A for the basis selected 
here would have been 32, for a grand total of 611; thus this basis, shown in 11.3A, would 
have been better for the problem of Sec. 10 than that of 10.2A. 

We may assume that the multiplications required to get the inverse from this form 
is the same as for the problem in Sec. 10.3. Thus 


11.38. Number of multiplications required for standard partitioning is 


577 + 1079 = |1656 


11.4. Standard inversion is still, of course, [4096 , 








22 J. PAUL ROTH [Vol. XVII, No. 1 


bo 
ie) 


4 § 6 7 8 9 1011 12 13 14 15 16 


ie =] x | 
| 


a  -! 








ph 

* 

A 
4 
A 








16 x | x x x | 





Fic. 10. JL’ for basis chosen. 


Appendix. One very important advantage of tearing over standard methods of 
solving networks is the fact that the solution developed for a component of the torn 
system is in fact a perfectly complete solution for this subsystem in itself. Thus regardless 
of how the subsystems are to be connected together, their individual solutions are 
completely useable. This is one advantage which solving a network by partitioning 
Y’ (or L’) for a given fixed system cannot enjoy. For the partitioning is performed on 
the Y’ for a given fired system; in the Y’ the interactions both through the topological 
connections and through the ‘electromagnetic couplings’? may be rather thoroughly 
scrambled: partitioning Y’ cannot sort these effects out. 

Another way of contrasting tearing with standard methods is this: A network problem 
consists of a pair (L, K), an ohmic tensor L is a complex K compatible with L. In tearing 
one seeks another network problem (L*, K*) in which the tensor L* (or Y’) splits in a 
convenient manner. From the solution of this network problem one then deduces the 
solution for the original problem. In solving the network problem (L, K) by standard 
methods, one operates only on L’ or Y’ for this given network problem; for instance, one 
partitions Y’ in the most expedient manner. Both of these advantages may be exceedingly 
important in the solution of very large physical problems. 
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The examples of this paper exhibit definite numerical advantages in tearing but, 
more importantly, they carry forward and improve on the mechanics of Kron’s method 


of tearing. 
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THERMAL INSTABILITY OF VISCOUS FLUIDS* 


BY 
CHIA-SHUN YIH 


University of Michigan 


I. INTRODUCTION 


The stability of a viscous fluid in an insulated vertical tube or between insulated 
vertical planes when a negative temperature gradient is maintained in the upward 
direction depends on the magnitude of this gradient, the gravity, the geometry of the 
solid boundary, the properties of the fluid, and the wave length of the disturbance. 
The purpose of this paper is to present the relationships between these variables for 
neutral stability, and the results concerning the effect of rotation on stability. 

The problem of stability of two superposed fluids in a cylindrical tube with surface 
tension at the common interface was solved by Maxwell in [1], without consideration of 
viscosity or wall effects. The stability of a layer of viscous fluid between two infinite 
horizontal planes when heated from below was investigated by Rayleigh [2], Jeffrey 
[3], Low [4], and Pellew and Southwell [5]. The problem under consideration resembles 
Maxwell’s in geometry and Rayleigh-Jeffrey’s in the means of producing instability 
as well as in the nature of the physical process through which stability is maintained 
under certain conditions. This problem has already been considered by Hales [11], 
Taylor [12], and Ostrach [13], but only incompletely. The contributions of these authors 
will be referred to at the appropriate places in this paper. 

As will be seen, the curves for neutral stability differ from those ordinarily obtained 
in investigations of hydrodynamic stability in that the critical Rayleigh number occurs 
at zero wave number, and hence that these curves have no lower branch. 

The “principle of exchange of stabilities’, which is assumed in many investigations 
of hydrodynamic stability but proved only in one instance (Pellew and Southwell, [5]), 
has been shown to be valid without general rotation. With the presence of general rotation, 
this principle is valid under certain restrictions. The investigation of the effect of rotation 
lends some support to the belief that rotation has no effect on the onset of instability. 


II. Srapinity oF Ftuip BETWEEN PLANE WALLS 


1. Formulation of the problem. If a layer of heat-conducting viscous fluid between 
two vertical planes is heated from below, free convection will occur only if the (negative) 
temperature gradient in the vertical direction is sufficiently great in magnitude. The 
stability of such a fluid layer is discussed in this section. 

Using Cartesian coordinates (x, , 2 , 23), with z; measured in the upward vertical 
direction, one can write the equations of motion as 


0 ou; 


ou, Ou;) _ dia, os aa 2 ( a) “= 
(% + u; 2) = (0,0, —g)p az, + "az, \Paz,)? = 12 3) (I) 


in which the summation convention has been used, p is the density, 7 is the time, g is 
the gravitational acceleration, p is the pressure, and » is the kinematic viscosity, which 
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is assumed to be constant. The velocity component in the 7th direction is denoted by 
u,; , and the three quantities in the parenthesis on the right-hand side of Eq. (1) are 
associated with indices 1, 2, and 3, respectively. The equation of continuity is 


0 Ou; 
ce U; =f —- = 0), (2) 
OT Ox; Ox; 


and the diffusion equation is 


oT a) 0 ( ar) 
<= + 4, —) = c——(p—}, (3) 
Az Fi 5e,) = "az, \? az, 
in which T is the absolute temperature and «x the thermal diffusivity, assumed to be 
constant. For moderate temperature differences the equation of state can be approxi- 


mated by 
p = poll — a(T — T,)], (4) 


in which 7, and pp are the temperature and density, respectively, of the fluid at a point 
chosen to be the origin, and a is the thermal expansivity. 

If convection is present, the temperature, pressure, and density will differ from 
their mean values. One can write 


T=T.+7", 
P=PDant?D’; (5) 
P = Pn t+ p’, 


in which the subscript m indicates primary quantities and the primes indicate the per- 
turbation quantities. For the primary temperature distribution 


re = T > + B23 ] (6) 
one has, from Eq. (4), 
Pm = poll — aB2s), (7) 


so that the hydrostatic pressure is given by 


_ = —gp(l — aBz3). (8) 

The thermal expansivity for water under normal conditions is of the order of 0.0001/°F, 
that for air is of the order of 0.002/°F. Thus if the maximum temperature difference is 
not excessive, the actual change in density is small. It will be assumed here once and for 
all that the only effect of density change is on the body force per unit volume due to 
gravity—the effect on the inertia or specific heat capacity being neglected. Although 
the change in specific weight is small because a is small, it must in no circumstances be 
neglected, because this change is the motivating force of any convection, and the sole 
cause of instability. 

The requirement of small change in density imposed a limitation to the magnitude 
of the maximum value of z; or of x; . If later one does not hesitate to speak of “zero 
wave number,” which corresponds to infinite wave length in the 2;-direction, it will be 
with the understanding that that term is only a convenient expression for “long wave 
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lengths.” As will be seen, the rate of change of the Rayleigh number with wave number 
at neutral stability is small when the latter is small, so that the result for ‘‘zero wave 
number” applies rather accurately for long wave lengths. 

Under the assumptions made on the effects of the density change, neglecting quadratic 
terms of the perturbation quantities, and remembering that the undisturbed state is 
one of dynamic and thermal equilibrium, one can write Eqs. (1) to (3) as 


Ou, Ee ioe 1 Op C 
ar = (0, 0, gal ) Po ax, + vAu,; , (9) 
Op’ Ou, _ 
Or aB polls + Po Ox, =? 0, (10) 
0 : 
(2 — ca) ‘= — fu, , (11) 
OT 


in which A is the Laplacian operator in Cartesian coordinates. By virtue of Eq. (4), 
Iiq. (10) can be written as ; 
IT’ Ou; 
~of 2 -+ au) + —! = 9Q, 
OT Ox; 
which, because of the smallness of the thermal expansivity, becomes 
oui _ 9 (12) 
—_ = 
In other words, the effect of change of density on continuity can be neglected. Eliminating 
p’ from the first and third equations contained in Eq. (9), one has, 
0 {odu, aus) oT’ (2 au 
— — —] = —-geae — A\ — — =—}- (13 
Or (2 Ox, ge zy 7? OX3 Ox, (13) 


For further discussion u, will be assumed to be zero. This does not mean that the 
motion under consideration is necessarily two dimensional in the usual sense of the word, 
for u, and u,; may still depend on 2, . If u, is zero, however, the equation of continuity 
permits the use of Lagrange’s stream function: 


“= U: 
, OX; ‘ . Ox, 


Thus Eqs. (11) and (13) can be written as 


a oY f 

(2 ra)t — *s. : (15) 
0 oT’ —- 

(2 _ ra)ay = ga = (16) 


If the origin of the coordinates is taken midway between the plates and the half spacing 
is denoted by d, the boundary conditions are 


9 te) - 
y = 0, < ees 0, and —=0 at 2x, = +d. (17) 
Ox; Ox, 


- 
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For convective motion independent of the coordinate x, and periodic in the z;-direction, 
one tries solutions of the form 


Y = «f(z) cos az e"’, (18) 
T’ = B d@x) cos az e”’, (19) 
in which 
_({% 2% 2s ia 
(«, y,2) = (2,2, %), ' a? 


a is the dimensionless wave number, and o (equal to o, + %i0;) is the complex amplifica- 
tion factor. Equations (15) and (16) now become 


[co — (D® — a’)]o = —Df, (20) 
[e — Pr(D’ — a’)|(D° — a’)f = —R Pr Dé, (21) 
in which D denotes differentiation with respect to z, Pr is the Prandtl number v/x, and 
4 
KV 


is the Rayleigh number. The boundary conditions are 
f=0, Df =0, and D@é=0 at x= +1. (23) 
If u, is zero but the motion is assumed to be periodic both in z, and 2; , one tries 
solutions of the type 


¥ = «f(z) cos by cosaze”, (24) 
T’ = Bd@(x) cos by cosaze”. (25) 
Equations (15) and (16) now have the form 
[co — (D’? — BW? — a’) Jo = — Df, (26) 
[o — Pr(D®? — b’ — a’)|(D® — BW — a’)f = —RPr DO. (27) 


The boundary conditions are the same as those for the strictly two-dimensional case, 
but the differential equations now correspond to truly cellular convection of a viscous 
fluid. Since Eqs. (26) and (27) would be identical to Eqs. (20) and (21) if a® + 6b’ were 
replaced by a’, the stability or instability against cellular convection can be predicted 
from that against the formulation of vortex tubes. This result is similar to that of Squire 


[6]. 

2. Principle of exchange of stabilities. In investigations of hydrodynamic stability 
other than that of the Tollmien-Schlichting type, it has often been assumed that the 
imaginary part of the factor o is equal to zero as well as the real part at neutral stability. 
This is the so-called principle of exchange of stabilities. Only in the case of thermal 
instability of a viscous fluid between horizontal plates has it been rigorously proved 
(Pellew and Southwell, [5]). This principle will now be proved for the problem formulated 
in the last section, but without any assumption concerning boundary geometry. 
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If the velocity and the pressure are assumed to be periodic vertically, and if 
u, = Ue", p’ = Pe", T’ = 6c", (28) 


the proof can be achieved by the use of Green’s theorems: 


I m,F, dS = I a dV, (29) 
> 


[| pe ag = [/ (grad F)-(grad G) dV + Il FAG av, (30) 
s V y 


on 


in which m; are the direction cosines of the outwardly drawn normal to the surface S 
enclosing a cellular space, and n is the distance along this normal. Multiplying Eq. (9) 
by u* , summing over 1, and utilizing Eq. (12), one has 


- 1 fpf Pus) i ane 
Jy = gol — — Il a + fff U*AU, aV, (31) 
Vv Vv 
in which 


a I U.UtdV, H= If oUt aV. 
V Vv 


The second integral on the right-hand side of Eq. (31) is, by Green’s first theorem, equal to 
/ [ mpur as, 


which is zero because on the solid boundary U* is zero and the quantity PU* is periodic 
in z, , the volume V being a cellular space of the convection. The third integral on the 
right-hand side is, for the same reasons and by virtue of Green’s second theorem, 


-/ff (grad U,)(grad U*) dV = —J, (say). 
Thus Eq. (31) can be written as 


ado tvd, = gaH. (32) 


If now Eq. (11) is multiplied by 7’* and integrated, one has, by Green’s second theorem 
and because of periodicity and the insulation of the wall, 


ol, + «xl, = —BH* (33) 


in which 
Io = If 66* dV, I, = If | grad 6 |’ dV. 
iJ. A 


From Eqs. (32) and (33) it follows that 


—B(o*Jo + vJ;) = galol, + «l;,), 
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or 


o(galyo + BJo) + gaxl, + BvJ, = 0. (34) 


o(gal, — BJy) = 0. (35) 


From Eq. (35) one concludes that if o; is not zero 6 must be positive. If so, from Eq. 
(34) o, must be negative and the fluid stable—as is also to be expected from the physical 
point of view. Thus, under the restrictions of periodicity in the vertical direction and 
of infinitesimal disturbances, the principle of exchange of stabilities is valid for a viscous 
fluid contained in an insulated tube and heated from below. 

If there are two horizontal planes intersecting the tube just considered, which are 
kept at constant temperatures to create a temperature gradient 8, periodicity in the 
Z3-direction is no longer necessary for the proof, because the integral involving the 
pressure in Eq. (31) now vanishes on account of the vanishing of the velocity components 
on a solid boundary. Since the temperature fluctuations on the horizontal plates are 
zero, the proof presented in the last paragraph remains valid. 

3. Solution of the differential system governing stability. The differential system 
consisting of Eqs. (20), (21), and (23) will now be solved. Since the principle of exchange 
of stabilities is valid, for neutral stability one needs only to consider the system (with 


h = R6). 


(D? — a’)h = R Df, (36) 
(D? — a’)’f = Dh, (37) 
f=0 and Df=0 at z= +1, (38) 
Dh=0 at x= +1. (39) 

Differentiating Eq. (36) and substituting Eq. (37) into the result, one has 
(D? — a’)*f = R D’f, (40) 

with the boundary conditions 

f = 0, Df = 0, (D? —a@’)*f=0 at r= 21. (41) 


Equations (40) can be written in the form 
(L? — RL — Ra’)f = 0, (L= D?—- a). (42) 
The solution of the indicial equation 


m*> — Rm — Ra’ = 0 


Mm = (n oo 1\(E) - (43) 


in which 


,§ Seam’ + 
, = ~~ -- . (44) 


t 
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Only one of the two signs need be taken, and either one can be taken. If the three roots 
of m are denoted by 


mM, =w, — a’ (¢ = 1, 2, 3) (45) 


which may be complex, the solutions of Eq. (42) are exponential functions with exponents 
+w,x. The form of Eq. (40) permits one to resolve the problem into two parts—in the 
one f is even, and in the other f is odd. The solutions for even f are cosh w,;x and those 
for odd f are sinh w,2. The secular equation obtained from the boundary conditions and 
determining the relationship between a and R is, if f is even, 
| cosh os] sh ws | 
| sh w, cosh w. Cosh ws | 
», sinhw, w,sinhe, w; sinh ws; | = ( (46) 
|m? cosh w, m: coshw. m3 cosh as | 
which, though complex as it stands, is only one real equation, on account of either the 
conjugacy of the two complex roots of m, or the fact that the three roots are all not 
complex. The secular equation for odd f is similar to Eq. (46), the only difference being 
that the symbols cosh and sinh are exchanged. 
It can be readily shown that 
Ow; Om; 
—=—— ion _— 
da da 


(¢ = 1, 2, 3) 


contain the factor a, hence vanish as a vanishes. To find the minimum Rayleigh number, 
one differentiates Eq. (46) or the corresponding equation for odd f with respect to a 
and sets the derivative dR/da to zero. Thus, for the critical condition and as far as the 
first derivative, R can be considered as a constant and all differentiations with respect 
to a can be considered as ordinary differentiations. After differentiating Eq. (46) or the 
corresponding equation for odd f, one obtains three determinants having one row with 
its three members containing the factors 

OW, OW. Jw, Om, OAM. AMs 

da’ da’ da ” da’ da’ aa 
respectively. Since these vanish for a equal to zero, the minimum Rayleigh number 
corresponds to zero wave number in the direction of gravitation—the possibility of a 
maximum being ruled out by physical considerations. 

The task of finding the minimum Rayleigh number is now very much lightened. 

One may set a equal to zero forthwith and solve the differential system directly. The 
equation to be solved is now 


(D° — R D’)f = 0, (47) 


with boundary conditions 

f= 0, Df = 0, D'f=0 at r= +1. (48) 
For antisymmetric motion (even f), the solution can be shown by a direct calculation 
to be 


—tan R'* = tanh R'™, 
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with 
; <R* <x. 
Thus 
Riv* = 2.365, R,, = 31.29 (49) 


for even f. After the writer found this number, he was informed by Dr. G. K. Batchelor 
that Sir Geoffrey Taylor already possessed it in 1953, though he never published it. 
The equation ieading to this number was also found by Ostrach for different boundary 


conditions. 
For symmetric motion (odd f), the solution is 


1/4 
tan R'’* = tanh R'* 


with 


Thus, approximately, 


53 J my 9" 2 
pis = — = 3.927, R., = 237.6 (50) 


Ver 4 


which, it must be remembered, is based on half of the spacing of the plates. This Ray- 
leigh number was given specifically as the critical one for the stated problem by Ostrach 






3 | ——o 
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Fic. 1, Approximate and exact neutral-stability curves for antisymmetric convection between plane 


boundaries. 
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Fic. 2. Approximate and exact neutral-stability curves for symmetric convection between plane bound- 
aries. 


[13], who considered only symmetric motion. Since it is greater than the critical Rayleigh 
number 31.29 for even f, the stability is governed by the latter. 

For non-zero values of a, the differential system consisting of Eqs. (36) to (39) will 
be solved first by Chandrasekhar’s method [7]. The results are shown in Figs. 1 and 2 
by dotted lines, for anti-symmetric and symmetric motions, respectively. After the 
approximate solutions are available as guides, Eqs. (45) and (46) are solved by trial and 
error. The results are given in Table 1. The corresponding curves are drawn in Figs. 




















TABLE 1 
: d | 0 | 0.5 | 1 2.5 | 3 
Sewanee | 31.3* | 37 | 53.5 360 | 742 
R(oddf) | 237.6% | 251.5 | 288 670 | 





*Given before by Eqs. (49) and (50). 


1 and 2 for comparison with the approximate solutions. The power of Chandrasekhar’s 
method is amply demonstrated. 


III. Srapiuity or Fium 1n A TuBE 


1. Formulation of the problem. In this part the stability of a viscous fluid in a 
tube of radius b and heated below is considered. The effect of rotation, which for a 
horizontal layer of fluid has been discussed by Chandrasekhar [10], will be investigated. 

If (x, , Z, , Zs) are cylindrical coordinates, the linearized forms of the equations of 
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motion and the equation of heat diffusion are, under the assumptions stated at the 
beginning of Part II 


Ou, Ou, : ; Op’ Uy 2 OU» = 
pl 4 + Q- Sia 202u, — Tem ~ 3 = Si. (91) 


1 On Ox, Zi Xi OL, 
OU. , OUs 1 dp’ Us 2 Ou, ia 
p ( 2+ 2+ 20u,) = ——F + py au. — 3 +55), (52) 
OT OX XY, OX x X) OX» 
OUs _ OU: Op’ 1 es 
pol 5 2 — ) = —f + povAus + pogal’, (53) 
OT OX» OX3 


in which the primes on p and 7’ indicate perturbation quantities, the equation of state 
has been utilized, and p, has the same meaning as in Eq. (4). The linearized diffusion 
equation is 

oT’ oT” , - 
—— — i > fe, « cit” (54) 
OT Ox 


in which £ is the primary temperature gradient in the direction of the vertical. If now 


one sets 
0) . r3 ’ 
7 »t25, = (1, ¢, 2) 
b b 


bu, bu. bus, 
(mm = ¥  -e = (u, v, w) 


K K 


b° b*p! 
T’=pb0, r=a—, qu, 


x? Po K 


the linearized equations become, with B for 2b*/« 


Ou ou Og U 2 ov 
—-+ (4 ) = —_— Pi - 4-32) 5! 
ry I ag 2v as + Pr\ Au a ap)? (55) 
ov ,{ ov 1 dq . v 2 du sa 
— +B + Quy = —- 3 + Pr dv -—-3a+3-], (56) 
ot dg r 0g : r 0g 
Ow Ow oq 
ae lh A | 2 ——— € mo) 
a, +B ap a, + Pr Aw Pr ko, (57) 
00 00 
—+B—+uv= 58) 
ry E ap w = AO, (58) 


in which A is now dimensionless and £# is the Rayleigh number based on b. The equation 
of continuity is 
O(ru) av , Arw) 
an ee aw ie eam am @, (59) 
or 6) Oz 


After the case of no rotation has been discussed the effect of rotation will then be 
investigated only for the rotationally symmetric case. If (as it turns out to be the case) 
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rotation has no effect on the stability of the most unstable mode, the most important 
aspect of the problem is solved. For axisymmetric convection, Eqs. (55) to (59) become 


du _ op, _ _ 99 nf - 5) 
re 2Bv = - + Pr\ Au 3)» (60) 
Ov v 
ry + 2Bu = Pr( a - *) , (61) 
Ow 0g 
— = ——-+ Pr Aw — Pr RO, (62) 
ot Oz 
= +w = AO, (63) 
ot 
Aru) , A(rw) _ 
0 dz *, (64) 


in which 


“2 4 
a=D+ip+5, p=2, r= -8@. 


dz yr” VK 
Equation (64) permits the use of Stokes’ stream function y in terms of which the velocity 


components are 


“= i oy : w= -* Dy. (65) 


r Oz 


Eliminating g between Eqs. (60) and (62), one has 


@ _ ply — 1) ](2u_ 2) _ op ® 4 p 
| 2 Pr\ A ) i ) = 2B a + Prk De 


r/\\dz or 


which, by virtue of Eqs. (65), can be written as 


) ; 
E es Pia _ 1) l(a ie 4) ¥ _ op 4 PR DO (66) 
ot Tr r/Tr Oz 
since 
du dw_ (pl To ~ 4)¥. 
Oz or (p r D+ r oy ne (4 rj/r 
Equations (61) and (63) can be written as 
a = ne =O = = 
|2 Pi(a 3) h - dz r’ (67) 
A 
(2 _ ao = 1 py. (68) 
ot r 


If one tries a solution of the type 


¢ = W(r) cos az ec’, 


v = V(r) sin aze"', 


0 = Or) cos az ee’, 
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Eqs. (66) to (68) become 


(o — Pr L)LW = 2aBV + PrR DB, (69) 
(co — Pr L)V = 2aBY, (70) 
tie ~ £96 = * Drv, (71) 


in which 


L=D'+*D-4-a'=D+Dr-a’, 


(72) 
=D +i p-d =i pD-@ 
r r 
The boundary conditions are 
v=0, DY = finite, V = 0, De=0 at r=0, (73) 
Vv = 0, Dy¥ = 0, V =0, De=0 at r=1. (74) 


The differential system consisting of Eqs. (69) to (74) governs the stability against 
axisymmetric convection of the fluid column under rotation and heated from below. 
If there is no rotation, B can simply be set equal to zero. 

2. The principle of exchange of stabilities. By a procedure similar to that used in 
Part II, it can be demonstrated that for axisymmetric motion a time-independent 
solution exists if 

2,8 3,8 

aa > 4, orif a > 4 
which means that for any given wave number a time-independent solution exists for 
sufficiently small Rayleigh number and sufficiently weak rotation. Furthermore it has 
been demonstrated* that, for zero wave number and undamped motion, 
2 

o = —nB = — 2, (75) 
in which 27/n is the period of motion in the direction of ¢. Equation (75) states that 
the convection pattern progresses with angular speed 2. For axisymmetric motion, 
n = 0, and o, = O if the disturbances are undamped. Thus for neutral stability axi- 
symmetric motion is time-independent for zero wave number. 

3. Solution for the case of no rotation. Since for the case of no rotation a time- 
independent solution corresponding to neutral stability exists, the differential systems 
to be solved are (with h = —R@) 


L’>¥v = Dh, (76) 
Lh = ‘ Drv, (77) 


*The demonstration is omitted to save space. 
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and 

v= 0, DY = 0, Dh=0 at r=0,1. (78) 
The differential equation and boundary conditions for V can be simply satisfied by 
taking V to be zero. Differentiating Eq. (77), one has 


LDh = R(L + a’)¥. 


Thus Eq. (76) can be written as 
Liv = R(L + a’)¥, (79) 


with boundary conditions 
Vv = 0, D¥ = 0, L’v¥=0 at r=0,1. 
Since the indicial equation is exactly the same as that preceding Eq. (43), the funda- 
mental solutions of Eq. (79) are the Bessel functions 
J (twr), J ,(tw2r), J ,(iwsr), 


in which the w’s have the same values as in Part II (with the Rayleigh number based 
on b, of course). The secular equation obtained from the boundary conditions is 


| J , (tw) J (tw) J :(tws) 
| wi J o(tw;) waSo(tw2) wed o(tws) | = 0, (80) 
| m2J (iw) mJ (tw2) mJ ,(tws) 

since 


Ld (tw . , 1 : 
hie = iwd (iwr) — : J, (twr). 


The determinant equation in which the w’s contain FR and 4, is a single equation though 
in the form given it is complex, and is the solution of the problem. With precisely the 
same arguments as in Part II, one concludes that the critical Rayleigh number occurs 


at zero wave number. 
For zero wave number a, the secular equation can be shown by a direct calculation 


to be 
J (RY) J (Ri) + iJ (RYT (Ri) = 0, 


the first root of which is 
R™* = 4.611, or R = 452.1. (81) 


This number was first given by Hales [11], who considered only the axisymmetric case 
which, as will be seen, does not correspond to the true critical Rayleigh number. 

Although only for the axisymmetric case has it been proved that the most unstable 
condition is associated with a wave number of zero, this situation can be expected to 
hold even for the other modes of motion. If one investigates the stability at zero wave 
number for the other modes (not axisymmetric), the conclusion reached in Sec. 2 of 
Part II enables one to write Eqs. (62) and (63) as 


L.W = Re, (82) 
L,0 = W, (83) 
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in which 
aes ; 2 l n 
w= Wr) cos ng, 8 = Ar) cos ng, L= D+-D-s4 
r r 
since one may assume the z-gradient of g (the perturbation pressure term) to be zero. 
This assumption can be justified physically from the symmetry of the flow, which is 
upward as well as downward and is motivated as much by the buoyancy of the hotter 
fluid as by the negative buoyancy of the colder fluid. The boundary conditions are 


W = 0, Dé=0 at r=1. 
W and 6 non-singular at r= 0. 
From Eqs. (82) and (83) it follows that 
L.W = RW, 
the adequate solutions of which are 
JAR ‘r) and 7@J,(iR''r), 
a combination of which is to satisfy the boundary conditions 
W = 0, DLW =0 at r=1. 
The secular equation is 
JiR) JUR) + iJ2GR’) JAR) = 0, (84) 
in which 
Ji(x) = - T(x) + Jn-i(2), 
the primes denoting differentiation with respect to the entire argument. The equation 


preceeding Eq. (81) can be obtained from Eq. (84) by taking n to be zero. The first 
roots of Eq. (84) for integral values of n are given in Table 2. Higher roots for each n 


TABLE 2 


n 0 1 2 
R'4 4.611 2.871 4.259 
R 452.1 67.9 329 .1 


can be found, which undoubtedly correspond to neutral stability. In this case, as in tlie 
case of plane boundaries, what happens as these higher roots are crossed has not yet 
been rigorously investigated mathematically or understood physically. As far as the 
first roots go, the second mode (n = 1) is the most unstable, and axisymmetric dis- 
turbances are more stable not only than those of the second mode but also than those 
of the third mode (n = 2). The second mode is antisymmetric, and compared with the 
antisymmetric motion for the case of plane boundaries is more stable. This is not sur- 
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prising, because the hydraulic radius for the latter case is, for b = d, exactly twice the 
hydraulic radius for the circular tube. The number 67.9 was first given by Taylor [12] 
without proof. 

For non-zero wave numbers, the method of Chandrasekhar has been employed to 
solve the system consisting of Eqs. (76) to (78). The results are shown in Fig. 3 by the 
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Fra. 3. Approximate and exact neutral-stability curves for axisymmetric convection in a circular tube. 


dotted line. With the assistance of these results Eq. (80) is solved numerically. The 
results are given in Table 3 and are represented by a curve in Fig. 3. These values are 
in good agreement with those given by Hales [11] for a smaller range. Hales observed 
from his numerical data for the symmetric case but without analytic proof that the most 
unsteady mode corresponds to zero wave number. 


TABLE 3 








Testes 











ane: Wes we 


R | 452.1 | 528 | 759 | 1322 








4. Axisymmetric motion with imposed rotation. Although the time-independence 
at neutral stability has been demonstrated only under certain restrictions when a general 
rotation is present, indications for it are so strong that it can be assumed. Putting ¢ 
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equal to zero and eliminating V and @ among Eas. (69) to (71), one has, for axisymmetric 


motion, 
3 2 ‘2aB\ 
E + RIL + a’) + (728) hy al 


By a process entirely the same as that used in Part II, one can show that for any definite 
value of B/Pr the critical Rayleigh number occurs at zero wave number. Since the 
effect of rotation is manifested in the number (aB/Pr)* which contains a’, for zero 
wave number the general rotation has no influence whatever on the stability. The critical 
Rayleigh number is therefore still given by Eq. (81). Physically, this is understandable 
because at zero wave number there is no radial motion, which alone is inhibited by rota- 
tion. In fact, by similar reasoning one may conjecture that a vertical magnetic field 
has no influence whatever on the value of the critical Rayleigh number, since no mag- 
netic lines are crossed at zero wave number. 

For non-zero wave numbers Chandrasekhar’s method will be extended to be used 
for three equations instead of two. If V’ = (2aB/Pr)V, 6’ = R@, the equations to be 
solved are 


Ly =-=— | kg = Dé’ (85) 
2 
LV’ = ~ (722) ¥, (86) 
Pr 
L’6 = -r(2 Drv) (87) 
with boundary conditions 

V’=0 at r=0,1; (88) 

v= 0, D¥ =0 at r=0,1; (89) 

De’ =0 at r=0,1. (90) 


If the \’s are the zeros of the Bessel function J,(A), the forms 


0 = DF) AadolAat), (91) 
m=! 


V’ = DY Aldi (nr), (92) 


are appropriate from the standpoint of the boundary conditions on V’ and @’. If these 
expressions are substituted in Eq. (85) and the result is solved exactly together with 
the boundary conditions, one obtains, 


o 


v= Z. An| Bui (iar + C,rJ (iar) + a2, 1.00 | 


m=1 


(93) 
+> A | Bis. + C!rJ,(iar) — ne Ji(\7) 
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in which 


To(mG 


(¢ m <a = (Am ’ i) MM ? 


: td ol 1a) 
(B., , BL) = (C.,. , CL) = 
, j ‘ J , (va) : 


— oe — rane 
M* =z, +a’, G = ~J.(Ga) 2J (ia) + iaJ,(ia). 


Substituting Eq. (93) in (87) and expanding in a series of Jo(A,7), one has 


M, _$f, [=@B.P. C.Qs _ Xs, | | wheat O20e re | 
R A, = > Aef Nom . Nom M? + An Nom + of \, (94) 





with 
iJ (\u)J (Ga) 
M . 





P, = 
1 


Q, = —2aP, +a / ir’ J (tar) J o(Xpr) dr. 


0 


Substituting Eq. (93) in (86) and expanding in a series of J,(A,,r), one obtains 
Pr \’ - —i,,B.P Con r; 
A ( ) A’ = i Bd mm i 
\ oaB = 2 | Nim + Nim + x, 


m= (95) 
-.  [ -r,.BeP.+CS. 1 
+ 2 As] Nin . | 





in which 


a il rJo(iar)T(dwr) dr, 


Nim = — 3d (Am) 2(An)- 


Taking m and n to be | in Eqs. (94) and (95) and demanding that A, and A{ not be 
both zero, one arrives at the condition 














xiP, = x2%; (96) 
in which 
ian M, 4. CQ: — a’B,P, be MD 
a R No M?’ 
__ C1 = BP) | 
” No M; F 
Ce an EP x? 

o, = —! . 1D yf 1 ; 

Ni + 





Pr\? | C{S,— BP, 1 
gap. ancl (7) 141 a. 5 a 
. . I, 2aB + N 11 M; 
Equation (96) represents the first approximation to the relationship between a, R, and 
(B/Pr)’ for neutral stability. Its graph for any fixed value of B/Pr intersects the neutral 
curve for no rotation at the R-axis, and can be expected to lie to the right of the latter 


curve for non-zero values of a. 
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IV. CONCLUSIONS 


From the foregoing it can be concluded that 

(1) With no rotation, convection at neutral stability is time-independent. With 
rotation the time-independence of axisymmetric motion for neutral stability can be 
proved only under certain conditions, but is valid if the wave number is zero. For other 
modes of motion in a rotating circular tube the undamped motion is time-independent 
relative to a frame of reference rotating with the tube, provided the wave number is zero. 

(2) With no rotation the critical Rayleigh numbers occur at zero wave number for 
plane boundaries and for axisymmetric motion, and can be expected to occur at zero 
wave number for other modes of motion in a circular tube. Antisymmetric modes are 
more unstable. For such modes the critical Rayleigh number is 31.3 for plane boundaries, 
and 67.9 for a circular boundary. (The priority of these numbers belongs to Taylor.) 

(3) Detailed relationships between the Rayleigh number and the wave number at 
neutral stability are given by Eq. (46), one similar to it, and Eq. (80), and by the graphs 
and tables. An approximate relationship connecting the Rayleigh number, rotation 
parameter, and wave number is given by Eq. (96) for neutral stability in the presence 
of rotation. Since the effect of rotation on the mean orientation of the fluid has been 
neglected in this paper, the results concerning the effects of rotation are only approximate. 
The results for zero rotation are rigorous. 
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THE CONSERVATION EQUATIONS FOR INDEPENDENT COEXISTENT 
CONTINUA AND FOR MULTICOMPONENT REACTING GAS MIXTURES* 


BY 
W. NACHBAR (Lockhead Aircraft Corporation, Missile Systems Division, Palo Alto, Cal. 
AND 
F. WILLIAMS anp S. 8. PENNER (Caliyornia Institute of Technology) 


Summary. The equations for conservation of mass, momentum, and energy are 
derived for a set of independent, coexistent continua obeying the laws of dynamics and 
thermodynamics. The idea of a control volume and a control surface for each continuum 
is used in the analysis. The derived results are practically identical with relations obtained 
previously by Th. von Karman. 

A direct comparison is conducted between the continuum theory results and those 
obtained from kinetic theory by assuming that, for each of the species, the kinetic 
theory definitions apply. It is found that the new terms appearing in the conservation 
equations derived from continuum theory are precisely those which are required to 
make these equations identical with the results obtained from the kinetic theory of 
multicomponent, reacting gas mixtures. However, the continuum theory forms of the 
equations are not useful because they require knowledge of the transport properties 
for individual species in the mixture. 

I. Introduction. The equations for conservation of mass, momentum, and energy 
for a one-component continuum are well known and are derived in standard treatises 
on fluid mechanics [1-3]. On the other hand, the conservation equations for reacting, 


t 


multicomponent, gas mixtures are generally obtained as the equations of change for 


the summational invariants arising in the solution of the Boltzmann equation [4, 5]. 
One of several exceptions to the last statement is the analysis of von Karman [6] whose 
results are quoted in a recently published book [7, 8]. Since von Karman’s method of 
analysis has not been deseribed in detail, and since his results seem to differ from the 
classical relations through the occurrence of higher order terms in the diffusion velocities, 
it appeared worthwhile to re-examine this problem with some care. 

The objective of our investigation is the derivation of the conservation laws for multi- 
component, reacting, gas mixtures. To this end we invent a physical model consistent 
with continuum theory. Our model involves the idea of a multicomponent continuum 
composed of coexistent continua, each obeying the laws of dynamics and thermody- 
namics, a notion which was first introduced by Stefan in 18717. For an n-component 
gas mixture we presume the existence of n distinct continua within any arbitrary volume, 
continuum K corresponding to the chemical species K. We shall use the terms continuum 
K, species K, and component K interchangeably, it being understood that each of these 
phrases refers to continuum K of the coexistent continua as long as we are following 


*Received December 13, 1957. This research was supported by the United States Air Force through 
the Air Force Office of Scientific Research of the Air Research and Development Command under 
Contract AF 18(603)-146. 

tWe arrive at the model of simultaneous coexistent continua as the logical transcription to continuum 
theory of the fact that the entire volume is accessible to all of the different molecules in a gas mixture. 
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the derivation of conservation laws from continuum theory. It is apparent that each 
space point in the multicomponent continuum has n velocities vs(K = 1, 2, --- , n), 
one velocity for each of the coexistent continua. 

In Sec. II we present relevant definitions and basic mathematical relations, which 
will be used in subsequent sections. In Sees. III to V we treat, respectively, the equations 
for conservation of mass, momentum, and energy. Our results are considered critically 
in Sec. VI and are compared with the relations obtained from the kinetic theory of 
non-uniform gas mixtures. 

II. Definitions and basic mathematical relations. The multicomponent continuum 
is considered to be defined in regions of space, every point in a region being an interior 
point of the region. All properties of the mn continua, including the velocities 
v\(K = 1,2, ---, 7), are assumed to be described by functions continuously differentiable 
in all variables within the region. This statement will be said to define ‘‘continuous flow”’ 
for the multicomponent continuum. 

The conservation equations for continuous flow of species K will be derived by 
using the idea of a control volume r*() enclosed by its control surface o“(¢), and lying 
wholly within a region occupied by the continuum; here ‘‘?’’ denotes the time. The 
notation of cartesian tensors will be used*. Let 2,(¢ = 1, 2, 3) denote the cartesian 
coordinates of a point in space. In cartesian tensor notation, the divergence theorem 
for any scalar function belonging to the Kth continuum a*(z; , t), becomes 


[ an’ dco = [pak dr, (1) 


“¢ 


where n*; denotes the outward normal to the surface o” and a; represents the gradient 


of the scalar a*. For any vector function belonging to the Kth continuum, u*(x, , 0), 
we have 
K_K K . 
| ujn;dc = | u;.; dr, (2) 
Jak Jk 
with u, denoting the divergence of the vector u* . 
Consider that some property of the Kth continuum has a density per unit volume 
K 1 


equal to a” (z; , t), and let A*(¢) be the amount of this property contained within the 


: 


control volume 7*. Thus 


K K ry) 
A th = | a (x, , t) dr. (3) 
. y, K k ° . . : y K: 
For example, if a = p = the density of mass of species K, then A” is the total mass 
of species K contained within r*. The property a* has a density per unit mass of mixture 
Ke 
equal to 6” (z; , t) where 
K 5K ean 
a = pb ’ (4) 
with p = > x", p* representing the density of mass for the fluid mixture. 


ry ° ° K ° . ° > K 

Che derivative dA*/dt is defined to mean the time rate of change of A” as the volume 
r* and its surface o* move with the flow of species K. Consider that Eq. (3) holds at a 
time f, ; at time ¢, + At the particles in 7“ at x, will have been displaced to new positions 


*Repeated subscript indices imply summation over all allowed values of the indices. 
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. . ‘ ° Kr 1 K 
x! and will be contained within some new volume r°(t, + At) = 7’” enclosed by a surface 
K F (Ks /K K ew 7 K K orp C 
ao (ty + Al) a’; in general 7’* and o”” are different from r* and o”. Therefore 


A*(to + At) = [ a (xi , to + Al) dr 


vr 


d . . | 
( aA" ) _* lim +2. if. a (zi, t + At) dr — , aX(x; , to) ar}. (5) 


\ 


and 


It is demonstrated in the appendix that Eq. (5) is equivalent to the relation 


a [ EE + (a*v’) Ja (6) 
dt — J-k ot hades = 
Hence, using the divergence theorem given in Eq. (2), it is fouud that 
1A*™ Aa* 
: - = [ * fer dr + I, a*y* a; do. (7) 
( dr 


Equation (7) expresses the idea that the time rate of change of A“ in a flow, for an arbi- 
trary volume r* bounded by a surface o“, is equal to the stationary rate of change of A* 
in the interior of r* plus the rate of change of A* due to the movement of r“ and o*. 
may be rewritten in the equivalent form 


d ; al og* bane 
dt (f ps" ar) ” I, E2 (o6") | dr. (8) 


The overall’ transport equation for the multicomponent continuum is then obtained by 


Mquation (6 


summing over components, a procedure which is in accord with the idea of inde >pende nt 
coexistent continua. We choose at the arbitrary time ¢ all of the control volumes r* to 
be coexistent, i.e., 7* = 7 for all K. We will henceforth refer to a volume r thus defined 
as being “of the multicomponent continuum” at time ¢. After summation, Iq. (8) now 


pecome 


~ fd 4 x F Gd , . ee. \ 
> (= I. p> ar) / | ( . ob") 4 x (08). | dr. (9) 


III. Continuity equations. We denote by w* the net production of mass of species 
K per unit volume per unit time. Since mass is not created or destroyed by chemical 


r 


‘tions, but only converted from one species to another, it follows that 
- 
z w = 0. (10) 
K 
Phe continuity of the mass of species K in an arbitrary volume r* is therefore expressed 
by the equation 
d 3 7K K 
- | pY° dr} = w drt, (11) 
dt \J_x Jk 
rK - . . . : rr K >K rK 
where Y“ is the mass or weight fraction of species K (i.e., p» = pY” and Y” equals the 


mass of species A in unit mass of mixture). Let 8“ = Y* in Eq. (8); then Eq. (11) becomes 


rs - > Ry 
| ES + (pY*v*) , - o*| dr = 0, 
K Cc 


vr 
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and, since 7* is arbitrary, we have 

* 7Ry 

K a(pY™) KK 
w = on (p “oF . (12) 
ot 

Now let v* , the flow velocity for species K, be represented as 

K 7K « 

v; =v+ ) es (13) 
where 


v= p ry. (14) 


The summation in Eq. (14) is extended over all n distinct chemical components. Thus 
vis the mass-weighted average velocity of the fluid mixture, and V“ is said to denote 


the diffusion velocity of species K. Since 
> Y* =1, (15) 
K 


it follows from Eqs. (13) and (14) that 
2 F"F7> = @. (16) 
AK 


Introducing Eq. (13) into Eq. (12) leads to the following equation for continuity of 


species K: 


: D oK\ - ae ae ee 
w = 2 (oY) + pl. + (oY V9), (17) 
where 
D 0 
Di ‘ = at? + i): U8) 


is the Euler total time derivative following the mass-weighted average motion of the 
multicomponent continuum. Summing Eq. (17) over all distinct components, in view 
of Eqs. (10), (15) and (16), leads to the overall continuity equation 


Dp 


- ys = 0. (19) 
Di + pv; 


Equation (19) is evidently also the correct form of the continuity equation for a one- 
component system. 

We may now transform Eq. (9) by using Eqs. (18) and (19) to obtain a form which 
is useful for the derivation of the differential equations expressing conservation of 
momentum and energy, viz., 


ri [ ps" dr) = [| > PE + © rev y..| dr, (20) 
K Qt J, 7X Jr K 


rT 


where 


B= DB". 
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IV. Momentum equations. For an arbitrary volume 7 of the multicomponent 
continuum, the total rate of change of linear momentum in the jth coordinate direction 
must equal the sum of the following: (a) the surface integral of the stress vector > x on, 
where of equals* the component in the direction x; of the stress vector acting on that 
face of an elemental parallelepiped of species K which has outward normal in the direction 
r, ; (b) the volume integral of the total vector body force Zs pf* acting on unit volume 
of mixture, where {* is the vector body force per unit mass of species K; and (c) the 
volume integral of the total rate of generation of momentum in unit volume through 
production of species. Let the rate of generation of momentum in unit volume for species 
K be w*m* , where m* is the average momentum of the generated mass of species K 
per unit mass of species AK. We postulate that, overall, linear momentum is neither 
created nor destroyed by chemical reactions; the consequent conservation principle 
states that the total rate of generation of linear momentum per unit volume by chemical 


production of species is zero: 
> w*m* = 0. (21) 
K 
The total rate of change of linear momentum is then expressed mathematically by 
] sia : ai 
Zz (: [ pY*v' dr) = [ > okn; do + [ > e*f% dr. (22) 
K dt J,x rk ear v¢ K dr K 
Using the divergence theorem, Eq. (1), and the transport relation given in Eq. (20), with 
B* = y*y»* 
B= 0%, 
Iq. (2) becomes 
P De’ ee - CoK 
| | » +(p >, Y*V‘4I A dr = [ > (oh « + p*f5) dr. (23) 
K Jr K 


We now define o? , the diffusion stress tensor, as 


op = —p D Y*V'VS (24) 


K 


and f; , the vector body force per unit mass of mixture, as 
f= Dfiy". (25) 
K 


Since 7 is arbitrary, Eq. (23) then leads to an expression for overall conservation of 


7 


momentum, namely, 


Dv; Ov; om c . 
p= = p— 4+ p= Dok +07 + of; (26) 
Dt ot K 
If we define o;; as 
¢;; = _ op tou, (27) 
K 


*The species vectors o*, n; and pf* represent the sums of all forces which act upon species K and 
which move with the velocity of species K in the mixture. These definitions are used in Sec. VI to identify 


our results with the results obtained from kinetic theory. 
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then Eq. (26) is reduced to the well known form of the momentum equation for one- 
component systems. Therefore, o,; is the stress tensor and pf; is the body force acting 
on an elemental parallelepiped which is moving with the mass-weighted average velocity 
v’ . Furthermore, we can express o;; as a sum of partial stress tensors o#;/* 


K 
= >» et", (28) 
K 
where, from Eqs. (24) and (27) 
<i rKypK 17K 
of” =0;; — pY V;V;. 


Each stress tensor can then always be expressed as the sum of a mean pressure tensor, 
a viscous stress tensor, and a viscous diffusion stress tensor; thus, 





K K V,K D,K F " coat »K 
ot * = —p*5,; + ri" + 777", where p< = —4o%;", (29) 
o;; = —pi;; + iy + TH ; where p= —to;,;. (30) 
The total pressure p is the sum of the partial pressures p* for the different species, i.e., 
p= dp" (31) 
K 
and so, in view of Eq. (28), it follows now that 
a (32) 
K 
t= >» =. (33) 
K 
The equations of von Karman [6, 7] are obtained by using Eq. (30) in Eq. (26), viz., 
Dv} Vv D 
p Dt = aes Er + (7,5 + Tii),é + pf; ° (34) 


V. Energy equation. For an arbitrary volume 7 of the multicomponent continuum, 
the first law of thermodynamics states that: 
tate of increase of (internal plus kinetic energy) = rate at which work is done 
on 7 by (body forces plus surface stresses) + rate of inward transport of heat by 
thermal conduction through the surface o enclosing r + rate of generation of 
energy through production of species within r + rate at which work is done on 
material produced within r. 


Let u* denote the absolute internal energy of species K per unit mass of species K and 
let u denote the absolute internal energy per unit mass of mixture. Then 


“= p i (35) 
K 


The kinetic energy of species K per unit mass of species K is 4v‘v4 . The total rate at 
which work is done on 7 by surface stresses and body forces is represented as the super- 
position of the rate of work done on the individual continua by their own surface stresses 
and body forces. 

For the mass w* of species K, which is generated by chemical reaction in unit volume 
per unit time, the sum of (a) the internal and kinetic energy carried by this mass, and 
(b) the work done on this mass in unit time, is w*(n* + im*m*), where n* is the average 
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specific enthalpy of generated mass of species K*. We postulate that, overall, energy 
is redistributed among various states but is neither created nor destroyed by chemical 
reaction; the consequent conservation principle states that the total rate of generation 
of (absolute enthalpy plus kinetic energy) per unit volume by chemical production of 
species is zero: 


> w*(nX + 4mm‘) = 0. (36) 


The analytical expression of the first law of thermodynamics, subject to the funda- 
mental postulates of independence and conservation, is therefore 


p | f pY*(uX + 4v%v') ar | 


k rX* rker 
(37) 
c K KK K 
= >» | ony; do + [ pf iv; ar | - i > an; do. 
K o t oe OUdi<kékd 

Here q; is the heat flux vector for species K, taken as positive for outwards heat transport. 
Equation (37) can be transformed by the use of Eq. (21) with 

Qk 7K, K 7K. K K 

BX = Y*u% + 4Y“v'v5 
an 


B=ut wi +4 > Y*ViVS. (38) 


K 
Since r is arbitrary, the following differential equation for overall conservation of energy 


then results: 


D - 
-(u + dy'v’) + 3 


PD“? bo (> YXVSV) + Io DO (YuVS + FVD 
K K 


(39) 
=p Dy fvit+ Div. — Lagi: 
K K K 
Let h*“ denote the absolute specific enthalpy of species K, which is defined as 
K 
Y*n* = Yuk + _ (40) 


The absolute specific enthalpy of the mixture is then h = }>x Y“h* = u + p/p. The 
total heat flux vector is assumed to be expressible in the form 


Da = -AT.,, (41) 
K 


where 7' is the temperature, and X is the thermal conductivity. Using the definitions 
given by Eqs. (29), (30), (40) and (41), Eq. (89) can be written in a desired form, viz., 


p a (u + 4uios) + (o DD Y*A*V')., = pfivs — (pri).s + [Cris + ravi). 
~ 


+ OT) +e DY VG A DAG + eV 42) 


+ Me DY VIVIVS)  — 0 iD YVIV). 


K 





*The quantity »“ should not be confused with the total (average) specific enthalpy of species K 


which we denote by A*, as in Eq. (40). 
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Eq. (42) is practically identical with the energy equation derived by von Karman [6, 7]. 
When the diffusion velocities vanish, Eq. (42) reduces to a well-known form of the energy 
equation for one-component systems. 

VI. Comparison between the conservation laws derived for independent coexistent 
continua and the kinetic theory results for multicomponent gas mixtures. In order 
to show that the model of independent, coexistent continua represents correctly a real 
mixture of gases composed of different chemical species, we must compare the results 
obtained from this model with those of the kinetic theory of non-uniform gas mixtures. 
Quantities such as the density p, the mass-weighted average velocity v/ , and body force 
f, , have obviously analogous meanings in both the kinetic theory and the coexistent 
continua model. On the other hand, the precise kinetic-theory meaning of terms such 
as the stress tensor o* , the absolute internal energy per unit mass u*, and the heat 
flux vector q*, is not immediately apparent. In view of the known success of continuum 
theory for one-component systems, we shall identify the continuum theory properties 
os, u“, and q* for species K with their kinetic theory counterparts. Our proof then 
involves a comparison between the conservation equations obtained from multicom- 
ponent continuum theory, replacing continuum properties for each species by their 
kinetic theory definitions, and the conservation equations obtained from the kinetic 
theory of non-uniform gas mixtures. 

A. Definitions of kinetic theory. Let c*’" be the velocity of a particular molecule, m, 
of species K, and let V/*'" be the velocity of this molecule in excess of the velocity v* , 
which is identified in kinetic theory as the mean velocity of all molecules of species K. 
Then 


K,m\ ryK im 


Km K rpK m rK r7K,m K / 

c;'" =v; + Vi =v+V;+ Vi", (;'") =v; and (Vj"'") = 0, 
where the angle brackets indicate an average over all molecules of species K taken with 
respect to a distribution function appropriate for the mixture. 

From kinetic theory, we have the following definitions’ for the properties of species 


K in the mixture 


of = —p (Viv) (43) 
“ue = vor + rr = wore + - (44) 
q* — p* (AV yrs a "tee | ede = rire re re § ores. (45) 


- os RK ° ° ° e 
In Eqs. (44) and (45), u*’”, the total internal energy per unit mass of a molecule of 


species K, is expressed as the sum of (3V/“'"V/*""), the peculiar translatory kinetic 


energy per unit mass, and (2 the contribution of additional internal energy terms 
(rotational, vibrational, etc. per unit mass. We have then defined u* (u : and 
-K ae 

t= 


The corresponding definitions in the kinetic theory for the properties of the gas 


mixture will be denoted by the superscript 7’; these are”: the mixture stress tensor o,} 


where 


1See, for example, Secs. 2.31, 2.4 and 2.45 in [4]. 


*See, for example, Sec. 2.5 in [4]. 
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05 = — Le (Viem + V(ViN™ + Vi) 
; (46) 


=— os KVeovs — > oXVEVE , 
K K 


° . ° T 
the internal energy per unit mass of mixture u’, where 


pu” = Do [20%(ViR™ + ViVi + V')) + ptm 
K fod 
(47) 
3 1 > p XV _ “oe ™) +} de KV"V 7K i+ ) » pi* 
K 


K 


and the heat flux vector for the mixture gj , where 
G5 = De oI VEE™ + VIVE + VA) HVE + VD) 
K 


DAVE MVE MVE) + AV VET) + VEE) ae 
+ AVEVIVE + GEE) + RV. 


Using Eqs. (43), (44) and (45) in Eg. ( (46) » (47) and (48), the following identities are 
obt ined between the properties o,7 , pu Y ond q’, of the gas mixture and the properties 
of , p‘u* and q* of the individual species: 


Doi — De Vivi, (49) 
K K 

pu" = p*u* ‘+i de” P ViVi a (50) 
y q* Po DaiV5 + dP K u™Ve +3 de” Pere . (51) 


In each of the above relations, the property for the mixture is equal to the sum, over 
all species, of the corresponding property for the components plus various diffusion 
terms. The diffusion terms arise because the reference coordinate system for species K 
is taken to move with velocity vf + V* , which is the mass-weighted average velocity for 
molecules in species K alone; the reference coordinate system for the mixture on the 
other hand is taken to move with velocity v/ , which is the mass-weighted average 
velocity for all molecules in the mixture. 

B. Comparison of conservation equations. The continuity equation for species K, 
as given by Eq. (17), is readily seen to be identical with the corresponding relation in the 
kinetic theory for multicomponent gas mixtures if w* is the net mass rate of production 
of species K per unit volume by chemical reaction*. Explicit evaluation of w* requires the 
introduction of the laws of chemical kinetics*. 

The expression for overall conservation of momentum, Eq. (26), is likewise identical 
with the corresponding relation in the kinetic theory’, as a comparison of Eqs. (27) and 


(49) shows that ¢;; = o/. 


8See, for example, Eq. (4) of Sec. 8.1 in [4]. For chemical reactions, the right hand side of this equation 
does not vanish but is equal to w*/m*, where m* is the mass of a molecule of species K. 

4See, for example, Chap. I in [7]. 

‘See, for example, Eq. (7) of Sec. 8.1 in [4]. 
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To demonstrate the equivalence of the energy conservation equations, we rewrite 
Eq. (39) in the form 


D ok - >*Ktrr rK 
pj (L(Y *u™ + AY*VGV4) + 4oio'] 


+(p DY Y*u*V").; = ofvi + op D Y*fiv™ 
K K 


(52) 
+ Deis + Di. — Lai 
K K K 
— (pv; Ps Pr aa ~ Mes, EVI: - 
rE 
teplacing >> x c* , Dox Y*u* and >>, q* by their kinetic theory equivalents, as given 
by Eqs. (49), (50), and (51), Eq. (52) becomes 
Df er , 7K ,K 17K T Ty K 
P Di [u + 30705] - pf iv; . p >» y fi) oo Sia + (0,304). ° (53) 
t K 


Multiplying the momentum conservation equations by v/ and contracting, the following 
scalar equation is found: 


1 D ~ 
5 PT, Vivi) = 081.0) + ofr; - (54) 


Using this relation and Eq. (19), Eq. (53) can be placed in a form which is identical 
with the usual form of the result obtained from kinetic theory", viz. 


7 (pu") + puvi.s =p DL YPVS — ahs + ori. (55) 

C. Concluding remarks. We have shown that, with proper identifications, the 
conservation equations for mass, momentum, and energy which are derived from the 
model of independent, coexistent continua (multicomponent continuum) are those 
obeyed by real gas mixtures. This conclusion is not surprising, in view of the fact that 
the results derived from this model, as well as those obtained from kinetic theory, do 
not depend on the forces operative in molecular collisions. The total mass, momentum, 
and energy are conserved in collisions; they are summational invariants. One would 
expect, in general, that the independent, coexistent continua model will give the correct 
conservation equations for summational invariants. 

However, apart from extending the one-component continuum results in a natural 
way to the flow of reacting mixtures, the multicomponent continuum model does not 
lead to any new results which are presently useful. In particular, to express conservation 
of energy for the mixture, Eq. (42) requires knowledge of species transport terms u*, 
a, and q’; in the mixture; these terms cannot be evaluated by kinetic theory methods. 
The mixture transport terms u’, o,7 and q7 can be evaluated, however, and therefore 
it is Eq. (55), or an equivalent form, which must be used to express energy conservation 


for the mixture. 





*See, for example, the equation preceding Eq. (8) of Sec. 8.1 in [4]. 
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APPENDIX 


Proof of Eq. (6). In order to prove Eq. (6), it is more convenient to work 
from Eq. (3) than from the limiting relation given in Eq. (5), and also to introduce the 
Lagrangian representation [9]. For any continuum K, let the three parameters a* identify 
the individual point particles of continuum K; for definiteness, suppose that a% are the 
spatial coordinates of the particles of continuum K at some fixed time t, . The spatial 
coordinates x; for any particle a* at any time ¢, t > t , are then assumed to be given 
by the functions* 2*#*(a% , t) which are taken to be single-valued and at least twice 
continuously differentiable with respect to each of their variables: 


i = r¥*(a‘5 ’ i), t > ty , (A-1) 


and aX = z** (a* , t)). The transformations are assumed to be one-to-one, so that the 
inverse transformations a* (zx, , t) also exist and are twice continuously differentiable. 
The flow velocities or “particle velocities” for continuum K, v#*(a* , t), are then defined 


as° 


K;__K ax*™ K 
vt (a;,)= = v(x; a, 
ot 
where the v*% are defined by the inverse transformation. Similarly, the Jacobian of 


. ‘ ; K K 
Eq. (A-1) is given as A*” or A’: 





A**(a* , t) 


ll 


K 
0a; 


det {ae} = A*(z, , 1). (A-2) 


If the integral of Eq. (3) is changed with the use of Eq. (A-1) to an integration at time 


t t, over the volume r* , then 


A*(t) = | _at(as , DA** drs. (A-3) 

The definition of the time derivative given in Eq. (5) is therefore equivalent to 

dA* (2a aK #K oa) . 

— = — ——] dro. A-4 

dt [, a te "eT a6 
But it is readily shown that’ 

BAK _ UE jek oe oF aek 
at Ox; sine 


and, therefore, Eq (A—4) may be written as 
dA* da** 
Ef (Oot sata ws 


Transformation of Eq. (A-5) to spatial coordinates leads to Eq. (6). 
*In this Appendix only, an asterisk on any function indicates that its variables are a; , t; functions 


without asterisks have the independent variables z; , t. 
See [5], Eq. (7.07). 








54 


1, 


moo to 


or 


~J 


W. NACHBAR, F. WILLIAMS, AND §&. 8. PENNER [Vol. XVII, No. 1 


REFERENCES 
H. W. Liepmann and A. E. Puckett, Introduction to aerodynamics of a compressible fluid, John Wiley 


and Sons, New York, 1947, chap. 7 
S. Goldstein, Modern developments in fluid dynamics, vol. I, Clarendon Press, Oxford, 1938, pp. 95-105 
L. Prandtl, Essentials of fluid dynamics, Blackie and Sons, Ltd., Glasgow, 1952, chap. II etc. 


. 8S. Chapman and T. G. Cowling, T’he mathematical theory of non-uniform gases, Cambridge University 


Press, Cambridge, 1953, chap. 3 


. J. O. Hirschfelder, C. F. Curtiss and R. B. Bird, Molecular theory of gases and liquids, John Wiley 


and Sons, New York, 1954, chap. 7 


}. Th. von Karman, Sorbonne Lectures, Paris, France, 1950-51 


S. S. Penner, /ntroduction to the study of chemical reactions in flow systems, Butterworths’ Scientific 


Publications, London, 1955, chap. 2 


. For related treatments, see the following: 


J. Stefan, Sitzber, Akad, Wiss. Wien 637, 63 (1871); 

C. Eckart, Phys. Rev. 58, 267 (1940); 

C. Truesdell, J. Rational Mech. Anal. 1, 125 (1952) 
R. Courant and K. O. Friedrichs, Supersonic flow and shock waves, Interscience Publishers, 
York, 1948, Sec. 7 


New 











A MODIFICATION OF PRAGER’S HARDENING RULE* 


BY 
HANS ZIEGLER 


Eidgendéssische Technische Hochschule, Ziirich 


1. Introduction. The response of a rigid-work-hardening material can be described 
by 

a) an initial yield condition, specifying the states of stress for which plastic flow 
first sets in, 

b) a flow rule, connecting the plastic strain increment with the stress and the stress 
increment, 

c) a hardening rule, specifying the modification of the yield condition in the course 
of plastic flow. 

It is customary to represent the yield condition as a surface in stress space, convex 
[1] and initially containing the origin. The current yield conditions for a metal are 
those of v. Mises [2] and of Tresca [3]. The flow rule generally accepted [4, 1] is also 
due to v. Mises [5]. It states that the strain increment vector lies in the exterior normal. 
of the yield surface at the stress point. As to the hardening rule, there are mainly two 
versions in use. The rule of isotropic work-hardening given by Hill and Hodge [6, 7] 
assumes that the yield surface expands during plastic flow, retaining its shape and 
situation with respect to the origin. Another rule, developed by Prager [8], assumes 
that the yield surface is rigid but undergoes a translation in the direction of the strain 
increment. 

The rule of isotropic work-hardening does not account for the Bauschinger effect 
observed in the materials in question. Prager’s hardening rule accounts for this effect. 
However, as Perrone and Hodge Jr. [9] have shown in special cases and Shield and the 
author [10] in a general investigation, Prager’s hardening rule is not invariant with 
respect to reductions in dimensions possible in almost any applications. In other words: 
if the yield surface in 9-space o;, moves in the direction of the exterior normal at the 
stress point, the two-dimensional yield locus, e.g., in plane stress ¢, , ¢, does not do so. 
In certain cases, e.g., if only o, and 7,, are different from zero, the Tresca yield locus 
in the plane oc, , 7., even deforms. 

It is clear that the physical consistency of Prager’s rule is not affected by the phe- 
nomena last mentioned. However, they complicate the application of the rule particularly 
in eases which otherwise would be simple enough to lend themselves to a complete 
treatment. On the other hand, the investigations of Shield and the author have shown 
that, under v. Mises’ yield condition at least, the yield surface, in all special cases, 
according to Prager’s rule moves in the direction of the radius connecting its center 
with the stress point. This suggests that a corresponding modification of Prager’s rule 
might simplify certain problems. In the following sections such a modification is formu- 
lated, investigated and compared with Prager’s rule. 

2. The modified hardening rule. Let us consider an element of a rigid-work-harden- 
ing solid, referred to an orthogonal coordinate system x; . The state of stress of this 
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element can be represented by a stress point P in a 9-space o,, with origin O. In this 
space, the initial yield surface is represented by an equation 


F(o;.) = k° = const. (2.1) 


In the following, for simplicity, attention will be confined to initially isotropic materials 
for which the form of the function F is invariant with respect to a rotation of the stress 
state. An initially anisotropic material can be treated in an analogous manner. 

The hardening rule suggested by Prager assumes that during plastic deformation 
the yield surface moves in a translation. After a certain amount of plastic flow, it is 
given by 

Foi, — ax) = k’, (2.2) 
where the tensor a,, represents the total translation. Because a,;, is not necessarily the 
isotropic tensor 6;, , where 6;, is the Kronecker delta, the material becomes anisotropic as 
a result of the hardening process. Accordingly, direction is important, and we shall fix 
the coordinate system x; with respect to the element, small deformations being assumed. 

In the space o;, , the a,, represent the radius vector of the point C’, which initially 
was at the origin and which in the following will be referred to as the center of the yield 
surface. Due to the flow rule of v. Mises, the plastic strain increment de,;, , considered 
as a vector in the same space, lies in the exterior normal of the surface (2.2) at P. Thus, 
it is represented by 

iio. >. (2.3) 

OC i: . 

The definition of a Prager-hardening material is completed by assuming that the surface 
(2.2) moves in the direction of de,, ; more explicitly 


day, = cdéx , (2.4) 


where c is a constant characterizing the material. 
Instead of (2.4), let us assume 
daz = (oi bed Q ix) du, du > 0, (2.5) 


t 


i.e. that the yield surface still moves in a translation, but in the direction of the vector 
CP connecting the center of the yield surface with the stress point (Fig. 1). This rule 


4] 
YU 


Fic. 1. Hardening rule and flow rule for a linear work-hardening solid. 


is a modification of Prager’s law, physically acceptable since both sides of (2.5) are 
tensors of the second order. 

The scalar du in (2.5) is determined by the condition that P remains on the yield 
surface in plastic flow. If in 9-space the summation convention is adopted, this condition is 
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(do; = da;,) oF = 0, (2.6) 
OC +n 
and from (2.5) follows at once 
_ __(0F/de,;) doi; (2.7) 





b= > . 
_- (ox, = Oy) OF /dox, 


Since (2.4) has been replaced by (2.5), the vector OC no longer represents the total 
strain. This is a serious drawback of our modification of Prager’s rule. On the other 
hand, we gain a substantial advantage: dd in (2.3), i.e. the magnitude of the strain 
increment remains free and can still be established as a suitable function of the stress, 
the stress increment and the stress history. In other words: it is possible to adjust the 
modified hardening rule (2.5) to any kind of hardening law in simple tension and compression. 

The simplest way to dispose of dX, i.e. to complete the flow rule, is to assume that 
the vector cde;, is the projection of do,, (and thus of da;,) on the exterior normal of the 
yield surface. This corresponds to the procedure familiar from Prager’s rule, and it will 
turn out later on that on this basis the results of either rule coincide in many cases. If 
dd is fixed in this way, the total strain is represented by the sum of the infinitesimal 
translations of the yield surface in the direction of its exterior normal at the stress point. 





Since 
oF 
(doi, =e de;,) — ae 0, (2.8) 
OC“ 
we obtain from (2.3) 
IF /de;;) do;; 
ani _ ieee (2.9) 


c (OF /de;,) (OF /dcx,) 


By means of the hardening rule (2.5), (2.7) and the flow rule (2.3), (2.9) the material 
is completely defined. It is easy to see that these rules are a generalization to complex 
states of stress of a linear hardening law in tension and compression, Fig. 2, which 
exhibits a Bauschinger effect. Moreover, by assuming that c is not a constant, but a 





Fig. 2. Linear work-hardening in simple tension and compression. 


suitable function of the stress history, e.g. of the distance OC or of the dissipation work, 
it is possible to adapt our rules to any material with a given non-linear hardening law 
in simple tension and compression. 

For most of the following considerations it will not be necessary to restrict ourselves 
to linear work-hardening materials, i.e. to constant values of c. 
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3. Some general properties. In an initially isotropic solid the yield function takes 


the form 


F(ox) ad G[I (ex, ly EAG:s), T;(o;)], (3.1) 
where 
I, . Cit; I, = BO 5505; ’ I; _ ATES (3.2) 


are the invariants of the stress tensor. If the initial yield is independent of the mean 


normal stress, 


Fox + Bou) = Flo), (3.3) 


where @ is an arbitrary scalar. When plastic flow has set in, the yield function becomes, 


on account of (2.2) and (3.1), 


Flex —_ Q;;.) == GL (cu — Qin); I (os spam Qin)» T3(o 4 — ax) ]. (3.4) 


From (3.3) follows that the values of (3.4) remain unchanged when o,, is replaced by 
ix + Bd; : tf yield is initially independent of the mean normal stress, it remains so. 
On account of (3.4), the flow rule (2.3) becomes 


(22 Ol, dG oly OG ah ) dx. (3.5) 


ol, OCs ol, OO + Ol; OC i 


de, — 
Carrying out the differentiation of the invariants, we obtain 


0G . 0G 0G acca 
dex, = | — be + (On — an) + = (6:5; — @:;)(0;, — a) | Ad. (3.6) 
ol, ol, “ 
Let us assume now that the physical coordinate system initially coincides with the 
principal axes of stress. Then we have first 


on = 0, («#k) and ay, = 0. (3.7) 


From (2.5) and from (3.6) follows 
da;, = 0, («{#k) and de, = 0, (t 5 k). (3.8) 


Since the material is isotropic at the outset, this means that the strain increment tensor 
de,, and the tensor da;, are coaxial with the stress tensor o,, . The relations (3.8) remain 
ralid if the second assumption (3.7) is replaced by the weaker assumption 


a, = 0, (¢ x k). (3.9) 


It follows that, if the principal axes of stress remain fixed in the element from the start, 
the strain increment tensor and thus the strain tensor remain coaxial with the stress tensor. 
If the principal axes of stress rotate, (3.8) holds only in a first step, provided the 
principal system of stress is used as the physical coordinate system. If (3.8) shall hold 
in a second step, the coordinate system must be rotated between the first step and the 
second one. This rotation, however, violates (3.9): due to the anisotropy caused by strain 
hardening, the strain increment tensor is in general not coaxial with the stress tensor. 
Many problems of practical importance can be treated in a stress space of less than 
9 dimensions. In certain cases, e.g., a 3-space defined by the principal stresses is useful. 
From our last result follows, however, that this 3-space is inadequate where the principal 


axes of stress are not fixed in the element. 
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The properties obtained in this section apply without exception also for a material 
obeying Prager’s rule [10]. 

4. Treatment in subspaces. On account of the symmetry of the stress and strain 
tensors, the problem may as well be treated in 6-space. It is convenient here and par- 
ticularly for the subsequent specializations to denote the physical coordinates by z, y, z, 
the stresses by o, , -*: , Ty:, °** , and the strains by ¢, , --+ , &, *** , Where the dots 
indicate cyclic permutations. 

In the new notations the yield condition (2.2) reads 





ie ‘ 
F(oe, id a, ’ phates ’ Ty ess Qys ’ es ’ Tey 7 Ary io © ‘) aid k ? (4.1) 
where 7,, , T-, , *** have to be considered as independent variables. The flow rule (2.3) 
becomes 
OF oF oF 
de, = — dy, -:+ , de,, = =— DM, >>> , de, = =—- DO, °*:, (4.2) 
do, OT es OV xz 


and the hardening rule (2.5) takes the form 
da, = (a, — a,) du, +--+ , da,, = (ty, — ay.) dp, +: | 
. (4.3) 
da,, = (ty — Oy) dp, -* J 
Treatment in 6-space, however, requires the elimination of the stress components 


Ty (= Tye), *** , Of the strain components e,, (= ¢€,.), *** , and of the displacements 
- , which, on account of (4.3), are equal to the displacements a,, , +> . On account 


Qua» ** 


of the symmetry of the stress tensor 
Plo, 5 *** 9 Fenn °** 9 Cen 9 °°) © Jn *** 9 Fan °°) (4.4) 
Thus, the yield surface in 6-space is given by 
fon — Ge °** 5 Tye — Aye» *°*) (4.5) 
’ 2 
= Flo, — Gey °** 5 Tee > Ges °** 3 Tey — Og J BH. 
From (4.2) and (4.5) we obtain 


de, = xi dy, -*:, dy, = 2 de,, A. dd, aie. (4.6) 


Oo; OT,: 
This is the well known result that the flow rule of v. Mises remains valid in 6-space, tf 
the state of strain is represented by the engineering components €, , *** , Yue 5 
Further, from (4.3) follows 


da, = (¢, — a,) du,-::, day, = (Ty: — Gy.) du, ++. (4.7) 


Thus, the hardening rule (2.5) applies without modification also in 6-space. 

In many practically important cases, some of the stress components are identically 
zero. Starting once more in 9-space, let us denote the stress components present by 
o!, , the zero ones by o/{ . The initial yield condition is then 


F(oi, , off = 0) = Hoi.) = Fk’. (4.8) 


If we are not interested in the strains ¢/{ corresponding to the zero stresses of , we may 
treat the problem in a subspace o/, . Here H(c/,) defines a new yield surface. 
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On account of the hardening rule (2.5) 
‘y = 0. (4.9) 


Qik 
Thus, the yield surface, after plastic flow has set in, is given by 
, 44 “= 0) = A(o!, — a},) = F’. (4.10) 


J , 
Fict, ™ Ges oa ais j 


It follows that also in the subspace a‘, the yield surface moves in a translation. 
On account of (4.10) 


0 ” 
de}, = ci dX = cf dx. (4.11) 
OC ix Oss 


Thus, the flow rule remains valid in any subspace. Of course, it supplies only the strain 

components ¢’, defined in this subspace, although the «¢/{ , too, may be different from zero. 
Finally, from (2.5) follows 

dai, = (ots me air) du, (4.12) 

i.e., also the hardening rule (2.5) remains valid in any subspace. This is another advantage 

in comparison with Prager’s rule which in most subspaces applies only in a modified 


form [10]. 
In the next sections, we shall discuss a few special states of stress particularly import- 


ant in applications. We shall restrict ourselves to materials where yield is independent of 
the mean normal stress. Here, as in (3.3), 
Ko — Gy HB, *** 5 Fen — yn 5 ***) = Nog — Gay *** y Ten — Me ***)s (4.18) 
5. Plane strain. Here, 
Ty: = 7,3 = O, e, = 0, (5.1) 
0. Thus, the yield function has the form 


(5.2) 


by definition. From (4.9) follows a,, = a,, = 

Pos — Gg 5 Gy — By 5 Gy ~ Oe y Tay — Mey) 
On account of (5.1) and the flow rule (4.6), ¢, must be absent from (5.2). Using (4.13), 
we finally obtain the yield condition 


Ro, — C8 5 Oy — 0% , Toy — Coy) = ke’, (5.3) 


where 


aj=a,—a, and alj=a,—a,. 


If the material obeys v. Mises’ yield condition, we have initially 


" - e x . 
(o, — o,) +47, = 3 7 ; (5.4) 


where a, is the yield limit in simple tension or compression. The yield surface therefore 
is the elliptic cylinder 
( , r\ 12 \2 4 2 — 
[((o, — at) — \o, — af) + 4(7., — a) = 3 70 (5.5) 
with its axis parallel to the line bisecting the angle of the first quadrant and with semi- 
axes (2/3)'a,, (1/3)'oo. 
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If Tresca’s yield condition holds, we have, instead of (5.4), 
(eo, — 0, +47, = a. (5.6) 
The yield surface is again an elliptic cylinder 


[(o, — al) — (0, — af)? + 4(7,, — @,,)° = a (5.7) 


with its axis in the same direction, but with semiaxes o)/2 °°, oo/2. 


Subcase a: If r,, = 0, (4.9) yields a,, = 0. The problem can be treated in a plane 
o, , ¢, , and the yield locus is the strip (Fig. 3) obtained by bisecting the cylinder (5.5) 
or (5.7) parallel to the plane a, , a, . 


\ p y 7" Gy (7resca) 


Fia. 3. Yield locus in plane strain with 72, = 0. 


Subcase b: If ¢, = 0, (4.13) can be used once more, and the problem can be treated 
in a plane og, , 7,, . The yield locus is the ellipse (Fig. 4) obtained by bisecting the cylinder 
(5.5) or (5.7) parallel to the plane o, , 7,, . Its semiaxes are 20,/3'", o,/3'”* for v. Mises’ 
yield condition and oo , oo/2 for the condition of Tresca. 





Fia. 4. Yield locus in plane strain with o, = 0. 


Since the length of the cylinder (5.5) or (5.7) is infinite, it does not matter whether 
its translation is given by the vector de in the direction CP or by its projection on the 
exterior normal. It follows that for a linear work-hardening law [c in (2.9) constant] 
the results obtained here coincide with those supplied by Prager’s hardening rule. 

6. Plane stress. Here, by definition, 


o, = Ty, = T,, = O. (6.1) 
From (4.9) follows a, = a,, = a, = 0. Thus, the yield surface has the form 


g(a, — Gy 5 Ty — Oy y Tey — Oy) = k’. (6.2) 
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62 
If the material obeys v. Mises’ yield condition, we have initially 
o, to, — 0,0, + 372, = (6.3) 


and thus 


9 


° + (¢, — a,)” — (0, — a,)(o, — a,) + 3(7., — a, = oH. (6.4) 


io. = & 
= ° . . . . . * o1/2 . . . . . ° 
The yield surface is an ellipsoid with semiaxis 2°’“a, in the direction of the line bisecting 
the first quadrant in the plane ¢, , o, , with semiaxis (2/3)'’*o, in the direction of the 
line bisecting the second quadrant and o,/3'”° in the direction of the axis r,, . 
If Tresca’s condition holds, we have initially the yield limits 


[(o, — o,)? + 472,]'” = oo (6.5) 


and 
4 | o. + 0, + [(o, — o,)* + 472,]'” | = 00. (6.6) 
Thus, the yield surface is an elliptic cylinder with the same orientation as v. Mises 
ellipsoid, closed by two elliptic cones. 

Subcase a: If r,, = 0, (4.9) supplies a,, = 0, and the problem can be treated in a 
plane co, , a, . If v. Mises’ yield condition (6.4) applies, the yield locus is the ellipse of 


Fig. 5 with the equation 
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Fic. 5. v. Mises yield locus in plane stress with rz, = 0. 


If the material obeys Tresca’s yield condition (6.5), (6.6), the yield locus is the hexagon 
of Fig. 6. 
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Fic. 6. Tresca yield locus in plane stress with 7, = 0. 
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Subcase b: If o, = 0, (4.9) yields a, = 0, and the problem can be treated in a plane 
o. , Tz, . The yield locus is the ellipse of Fig. 4 with semiaxes a) , o)/3'” in v. Mises’ 
case and a» , ¢)/2 in the case of Tresca. 

If P lies in a corner or vertex of the yield surface, the strain increment exhibits the 
same indeterminacy in direction as in a perfectly plastic material. In the lower right 
corner of the Tresca hexagon of Fig. 6, e.g., the vector de lies anywhere in the shaded 
region defined by the normals of the adjacent sides. However, P remains only in the 
corner or vertex if the yield surface moves with P, i.e. if the stress increment dé has the 
direction CP. 

7. Another special case. In certain cases, e.g. in a cylinder subjected to torsion and 


and tension, we have 
o, = o, = Tt, = 0. (7.1) 
l'rom (4.9) follows a, = a, = a,, = 0. Thus, the yield function has the form 
glo, — Q, y Ty: — Oy: » Ter — Azz) = k’. (7.2) 


If the material obeys v. Mises’ yield condition, we have initially 
a; + 3(ri, + Te) = 9 (7.3) 
and thus 


(o, — a,)” + 3(7,., — ay)” + 3(712 — a2)” = ~ (7.4) 


Che yield surface is an ellipsoid of revolution with semiaxes ao , 00/3”. 
If Tresea’s yield condition applies, the principal stresses are initially 


a= 0,  o,3 = 3[o, + fo; + A(t. + 72.)}"I. (7.5) 


The maximum shear stress is (¢. — o;)/2; hence, we have, instead of (7.4), 


9 


(o, — a)” + 4(7,, — ay)” + Ate — a2)” =o, (7.6) 
i.e. an ellipsoid of revolution with semiaxes oo , o0/2. 

Subcase «: If ¢, = 0, (4.9) supplies a, = 0. The problem can be treated in a plane 
Ty», Tez » Lhe yield locus is a circle of radius ¢)/3'”” in v. Mises’ case and o,/2 in the case 
of Tresca. 

Subcase b: If r,, = 0, (4.9) yields a,, = 0, and the problem can be treated in a plane 
o, , Tez - The yield locus is an ellipse as in Fig. 4 with semiaxes a , oo/3'”’ in v. Mises’ 
case and oy , ¢)/2 in the case of Tresca. It is clear that this result, apart from the difference 
in notation, is the one of Subcase 6b in Sec. 6. 

8. Discussion. By suitable linear transformations of the stresses (and of the a;,) 
the yield surfaces and yield loci obtained in Sees. 5 through 7 could be simplified. Ellip- 
soids can be transformed into spheres, ellipses into circles and so on. The hardening 
rule (2.5) is not affected by such a process. However, if the flow rule (2.3) is to apply in 
the new variables, it is necessary to transform the strains [10, 11] simultaneously. 

At the present stage of research, it seems hopeless to expect a decision by experiment 
between the hardening rules confronted here. Comparisons must be based, therefore, on 
purely theoretical reasoning. 

It has been pointed out in Sec. 5 that in plane strain the results obtained for a linear 
work-hardening material by the hardening rule (2.5) are the same as those supplied by 
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Prager’s rule. A comparison with [10] shows that the situation is the same in plane 
stress and in the stress state of Sec. 7, provided that the material obeys v. Mises’ yield 
condition. So far it does not matter which hardening rule is used. If, however, a material 
obeying Tresca’s yield condition is subjected to plane stress or to the state of stress of 
Sec. 7, the response depends on the hardening rule. Under Prager’s rule, the yield surface 
deforms in most cases. Under the rule (2.5), such deformations do not occur, since, 
according to a statement in Sec. 4, the hardening rule (2.5) applies without modification 
in any subspace. This is a definite advantage of the rule (2.5). 

Another point has been raised in Sec. 2: the possibility of adjusting the hardening 
rule (2.5) to arbitrary non-linear hardening laws in simple tension and compression. 
This advantage, however, is obtained at the expense of the geometrical interpretation 
of the total strain. If one is prepared to renounce this possibility, one clearly is in a 
position also to adjust Prager’s rule to non-linear laws. This point, therefore, is hardly 
of importance. 

A serious drawback of the hardening rule (2.5) lies in the indeterminacy of the strain 
increment in a corner or vertex of the yield surface. This problem arises already if the 
Tresca yield condition is applied to a perfectly plastic solid 

Let us assume that a perfectly plastic cylinder with axis z is subjected to simple 
tension. In plastic flow, the stress point has the position P in Fig. 7. The exterior normal 


Fic. 7. v. Mises and Tresca yield in simple tension. 


at P of v. Mises’ ellipse has two components with the ratio 2: (—1). Hence, the contrac- 
tion of the cylinder in the direction y is half the elongation in the direction x and thus 
the same as the contraction in the direction z, since the volume change is zero. In other 
words: if v. Mises’ yield condition applies, the deformation of the cylinder exhibits the 
same rotational symmetry with respect to the axis z as the state of stress. However, 
if the material obeys Tresca’s yield condition, the direction of de is free within the shaded 
angle formed by the normals in P of the sides 1 and 2 of the Tresca hexagon. In view 
of the incompressibility of the material this means that under Tresca’s yield condition 
the deformation may be arbitrarily unsymmetric. There is no reason why the cylinder 
should not contract, e.g., in the direction y alone. Experimental evidence seems to decide 
against Tresca’s yield condition and its associated flow rule. It is clear, however, that 
in experiments an exactly linear state of stress cannot be realized. In any event, if we 
accept the Tresca yield condition together with the flow rule of v. Mises, we have to 
admit the possibility of an arbitrarily asymmetric deformation in simple tension and 


compression. 
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The difficulty considered here is not present in a material work-hardening in accord- 
ance with Prager’s rule. Here, the strain increment in a corner of the Tresca hexagon 
is uniquely determined [10] by the stress increment, and the ratio de,/de, is 2: (—1) in 
simple tension, no matter which yield condition is used. However, a material hardening 
according to the rule (2.5) exhibits the same indeterminacy in de as a perfectly plastic 
material. In this respect the rule (2.5) is not quite satisfactory. 

We intend to compare the two rules in various examples. So far, it seems that the 
rule (2.5) is inferior to Prager’s hardening rule from a physical point of view, but easier 
to handle in certain classes of applications. 
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BOOK REVIEWS 


The physics of fluids. Edited by F. N. Frenkiel and a board of associate editors. Published 
by the American Institute of Physics. Vol. 1, No. 1 (Jan.-Feb. 1958). Subscription 
Price: $10.00 (U. S. and Canada), $11.00 (elsewhere), per annual volume of six issues. 


The new periodical will be devoted to original contributions to kinetic theory, statistical mechanics, 

structure and general physics of gases, liquids, and other fluids, including magneto-fluid dynamics, 
ionized fluid and plasma physics, rarefied gases and upper atmosphere phenomena, and liquid state 
physics and super-fluidity. Only original papers containing significant new results will be considered 
appropriate for the journal, but not purely expository papers. In addition to the full-length papers of 
the main body of the journal, there will be shorter Research Notes, as well as Comments and Errata. 
Whereas the new journal starts as a bimonthly, it is expected to become a monthly in the not too distant 
future. 
The first issue contains the following papers: “On the theory of the bubble chamber’’ by F. Seitz; 
“Wall effects in shock tube flow” by R. J. Eurich and D. B. Wheeler, Jr.; “Effect of radiation on shock 
wave behavior’ by R. E. Marshak; “Hydromagnetic stability of a conducting fluid in a circular mag- 
netic field’”’ by F. N. Edmonds, Jr.; ‘‘Statistical behavior of a reacting mixture in isotropic turbulence’”’ 
by S. Corrsin; “On the diffusion of a chemically reactive species in a laminar boundary layer flow’ 
by P. L. Chambré and J. D. Young; “Quantum statistics of interacting particles; general theory and 
some remarks on properties of an electron gas’’ by E. W. Montroll and J. C. Ward. 


Esercizi Di Meccanica Razionale. By Bruno Finzi and Paolo Udeschini. Libreria Editrice 
Politecnica Cesare Tamburini, Milano, 1958. viii + 544 pp. 


In the third edition of this collection of problems of rational mechanics for science and engineering 
students at the University of Milan, the authors have gathered a more complete set of examples taken 
from statics, kinematics and dynamics. The problems range in difficulty from the most elementary 
used in an introductory lecture to the fairly complicated encountered in the various fields of solid and 
fluid mechanics. 

All problems are entirely solved in literal form and are plane problems. No numerical examples 


are given. 
The book is a companion volume for the treatise on Rational Mechanics of Dr. Bruno Finzi, pub- 


lished by Zanichelli, Bologna, in 1950. 
Mario G. SALVADORI 


Integral equations. By F. G. Tricomi. Interscience Publishers, Ine., New York and 


London, 1957. viii + 238 pp. $7.00. 


In the opening lines of his preface, the author expresses himself as follows: ‘“‘One of the first subjects 
in mathematics to attract my attention was integral equations; yet this book appears after a score of 
others. Why is this? It is because the writing of a book on integral equations is a rather difficult task, 
a task for which many years of meditation are necessary. In fact, such a book must satisfy two require- 
ments which are not easily reconciled. In order to facilitate theoretical applications to existence proofs 
it must present the main results of the theory with adequate generality and in accordance with modern 
standards of mathematical rigor. On the other hand, it must not be written so abstractly as to repel the 
physicist, engineer, and technician who certainly need and deserve this mathematical tool.’’ It is the 
reviewer's firm belief that our gifted author has admirably succeeded in satisfying both requirements, 
and that this book is destined to take its rightful place among ‘‘the score of others’’ which the mathe- 
matical world has already been fortunate to receive from his pen. 

There are four long chapters, two brief appendices (one on algebraic systems of linear equations 
and the other on Hadamard’s determinant theorem) plus a useful collection of exercises, which is 


(Continued on p. 94) 
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THE EFFECT OF LEFT TURNS ON THE CAPACITY OF A 
TRAFFIC INTERSECTION* 


BY 
GORDON F. NEWELL** 


Brown University 


Abstract. A model is proposed for the traffic flow through a fixed cycle traffic 
signal on a narrow two lane highway. The average flows in the two opposing lanes are 
computed when the queues are arbitrarily long assuming that left turns are permitted 
and occur with fixed probabilities in the two lanes. It is found that the existence of left 
turns tends to favor short cycles of the light and, under certain conditions, the competi- 
tion between this effect and the obvious advantages of the long cycle light gives rise to 
an optimal cycle time at which the average flows have a maximum value. Some simple 
models for multilane highways are also considered and show optimal cycle times. 

It is possible to analyse in considerable detail the queuing problem for a fixed cycle 
traffic light on a single lane highway at which cars are not allowed to make left or right 
turns [1]. The corresponding problem for a highway wi'h two lanes of opposing traffic 
in which cars are allowed to make left and right turns (particularly the left turns) is 
much more difficult because cars that wish to turn left at a typical intersection with no 
special facilities for left turning cars interact with cars traveling in the opposite direction. 
The two queues of traffic in opposing lanes are therefore not statistically independent 
and one must consider the behavior of two queues simultaneously. 

One of the major complications in this problem arises from the fact that cars in 
lane | may turn left if there is a sufficiently large gap in the traffic of lane 2. To analyse 
this feature of the traffic flow, one must study rather detailed models which specify 
the size of the required gap and the statistical properties of the entering traffic. It is, 
however, relatively easy to determine from simple models the conditions under which 
the queues in both lanes will grow indefinitely. This defines the maximum traffic that 
the intersection can carry before both lanes have simultaneously reached capacity. If 
the traffic of lane 2 is below capacity, the capacity of lane 1 decreases with increasing 
flow in lane 2 and vice versa. The flows at which both lanes reach capacity is therefore 
also the minimum capacity of each lane as a function of the flow in the opposing lane. 

If the incoming traffic flows in both lanes equal or exceed the average rate at which 
cars can go through the intersection, the queues grow indefinitely and the behavior of 
cars at the intersection is independent of the queue lengths or the arrival times for the 
incoming flow. Furthermore, there are no gaps in the flow through which a car can turn 
left, except at the end of the green light period. The average rate at which cars can go 
through the intersection under such conditions therefore depends upon only a relatively 
few properties of the intersection. 

We consider here a model intended to simulate an intersection on a narrow two lane 
highway on which there is no room for passing at the intersection. Our purpose is to 


*Received Jan. 4, 1958. 
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find the flow through the intersection when capacity is reached or exceeded in both 
lanes. The model is described by the following properties: 

1. Each car in lane 1 will turn left with probability p and each ear in lane 2 will 
turn left with probability p’. 

2. During any cycle of the traffic light, cars will leave the intersection on a first 
come first served basis only at discrete and specified times t; < t, < --- < ¢, (the same 
for both lanes) with the most one car leaving the intersection in each lane at each time (; . 
For a less rigid interpolation of this, we could say that we observe the intersection only 
at discrete times ¢, , f. , --- chosen so that t; — t;_, represents some average time interval 
between cars leaving the intersection under the ideal situation in which all cars move 
forward. We then approximate the exact leaving times defined in some suitable way by 
a neighboring value of t; , either the preceding value of t; , the nearest t; , or the subsequent 
ralue of ¢; . It makes little difference how one actually does this. 

3. Forward moving cars or right turning cars in each lane will leave at consecutive 
times ¢; but a car that wishes to turn left must wait until a left turning car appears in the 
opposing lane or until the time ¢, , the end of the green cycle of the light. We neglect 
the possibility that a left turning car in lane 1 might be able to turn left if a car in lane 2 
turns right. This we can eliminate either by saying that there are no right turning cars 
or that this practically never happens even if there are right turning cars. Actually one 
could extend the model to include such possibilities but at the expense of complicating 
the algebra. Although the same techniques described below could still be applied to 
such an extension, the capacities would depend upon two more parameters, the prob- 
abilities of right turns in the two lanes. We prefer to sacrifice generality for simplicity, 
in this case. 

Although we expect this model to yield a qualitatively correct description of the 
traffic flow under the appropriate conditions, there are certain quantitative weaknesses 
in this model. Assumption ] is probably correct although there does not appear to be 
any experimental confirmation. Assumption 2 is a reasonable approximation for forward 
moving cars although it might be better to make the ¢; random variables. The main 
inaccuracy of this assumption, however, arises in conjunction with assumption 3 where 
we postulate that right turning cars leave at the same rate (measured in terms of the 
forward moving cars and left turning cars when opposed 


time intervals t; — ¢;_,) as 
we would expect that 


by left turning cars also leave at this same rate. In both cases, 
that the time intervals should be increased somewhat for the right and left turning cars. 
This model will therefore tend to underestimate slightly the influence of turns. The 
error is expected to, be relatively small, however, in most cases of interest because the 
main effect of the turning cars arises from the left turning cars waiting for several time 
intervals, a time which is probably large compared with the excess time required to 
execute the left turn once the opportunity arises. 

The mathematical problem of calculating the flow through the intersection can be 
solved by first observing that a “state” of the intersection can be defined at any time 
t; by specifying the desires of the first car in each of the two queues. From this one can 
valculate the probabilities of various states at any subsequent times. For all practical 
purposes, however, the right turn states and the forward moving states are equivalent 
in this model and will therefore be considered as the same states. States in which both 
cars move forward (or right) or both turn left are also equivalent since they both lead 
to the result that the two cars leave the intersection immediately and new cars appear 
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in both lanes for the next value of t; . We shall consider these states collectively as state 1. 
There are two other states. We denote by state 2 that state for which the car in lane 1 
wishes to turn left and the car in lane 2 wishes to go forward. State 3 will then be that 
state for which the car in lane 1 wishes to go forward but the car in lane 2 wishes to go 
left. State 2 leads to the result that the car in lane 1 remains at the intersection but the 
car in lane 2 leaves and a new car appears in its place at the next value of ¢; . 

We can define a 3 X 3 stochastic matrix P that gives the transition probabilities 
from any of the three states at time ¢; to any of the three states at time ¢;,, . This matrix is 


pp’ + (1 — p)l — p’) p’ Pp 
P = p(l — p’) (1 — p’) 0 . (1) 
(1 — p)p’ 0 (1 — p) 


The second column, for example, describes the transitions from the state with a left 
turning car in lane 1 and a forward moving car in lane 2. The car in lane 1 remains as a 
left turning car but a new car appears in lane 2. The final states are therefore states 
with two left turning cars (probability p’) or a left turning car in lane 1 and a forward 
moving car in lane 2 (probability 1 — p’). The third column describes a symmetric 
situation whereas column | gives the result of new cars appearing in both lanes. 

The first column of P also gives the state probabilities at time ¢, but, for uniformity 
of notation, it is convenient to represent the probabilities for the states at time ¢t, as 
given by Pe, with e, the vector (1, 0, 0). The state probabilities at time t; are given by 
the three components of the vector P’e, . 

For 7 < n — 1, the three state probabilities at time t; also describe the probabilities 
that cars leave the intersection at time ¢; in both lanes, in the second lane only, and in 
the first lane only. They also represent the expectations for these events. The expectations 
of the number of times these three events occur during the first 7 values of & are given 
by the components of the vector 


>> P’e, for j<n. 


Since at time ¢, , the cars in both lanes leave the intersection regardless of which way 
they go, the expectations of the total number of times per cycle that cars leave in both 


lanes, in the second lane only, and in the first lane only are given by the components 


of the vector 


n—1 n-1 


N e, + } P'e, = > P'e, : (2) 
k=1 k=0 
The sum of the first two components of N gives the average number of cars that 
leave in lane 2 per cycle, thus the capacity C, of lane 2 in cars per cycle (provided the 
queue of lane 1 never vanishes). Similarly, the sum of the first and third components 
of N gives the capacity C, of lane 1 in cars per cycle. The determination of C, and C, 
is thus reduced to evaluating the vector N. 
N can be evaluated most conveniently in terms of the eigenvectors and eigenvalues 
of the matrix P. If we let \, , A2 , and A; be the eigenvalues of P; f,; , f. and f; the right 
eigenvectors; and g; , g2 and g; the left eigenvectors defined by 


Pf, = ,f; and g;P = jg; , 
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then e, can be expanded in terms of the f; to give 
3 
oy = Dd flo » €1)/( Ge » te); 
k=1 


in which (g, , €,) denotes the usual real inner product of g, and e, . We then find 


3 n-1 
k=1 


7=0 


The sum over j is a geometric series, but since P is a stochastic matrix, one of the eigen- 
values (which we number \,) is 1 whereas | \, | and | A; | are less than 1. The summation 


over j gives 


; (9: ,@1) , (1 — A2)(qe , 1) (1 — A3)(gs , e 
N =n Gah + Owen! + Gada Bd @) 

The eigenvalues and eigenvectors of P can be found explicitly as: 
A, = l, A>, As = (1 — p)(l — p’) & [(1 — p)(L — p’)pp’)’”, (4) 
f Ip l—p)aA; +p’ — 1) | ; g p’(\; + p’ — 1)" (5) 

L'( — pA; t+p—1)- D(A; + p’ — 1)’ 
and so 

a lath ela f ° es ee PR * ne =a ™ 


Substitution of Eqs. (5) and (6) into Eq. (3) gives N in terms of p, p’ and n. The 
resulting explicit formula is rather awkward and it is easier to analyse the results in their 
present form. We consider below some special limiting cases. 

Case I. p = p’. The formulae for p = p’ are considerably simpler than the general 
case although in deriving the results it is somewhat easier to recalculate the eigenvalues 
and eigenvectors directly from P rather than using Eqs. (5) and (6) in which A; + p’ — I 


vanishes. The result is 


= oe. 
n (1 — p){l — (1 — p)"(1 — 2p)"] | 
a = 2 eee. eee. L <A a 
J 3 — 2p |! p| + n3 — 2p)? I (7) 
li—p fmigd 
From this we obtain 
Cag SS yg Oe 2 (8) 
(3 — 2p) p(3 — 2p) 


This formula is plotted in Fig. 1 for various values of n and p. In most practical cases 
we may consider n >> 1 (it is usually in the range 10 to 20). The important parameter in 
Eq. (8) for n > 1, however, is the quantity np which is nearly proportional to the total 
number of left turning cars per cycle. We can distinguish two cases. If n > 1 and np > 1 
(many left turning cars per cycle), then the exponential terms in Eq. (8) are negligible and 
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Fic. 1. The effect of left turning cars is shown by plotting Ci/n, the capacity C, with left turns relative 
to the capacity n with no left turns, as a function of n and the probability of left turns, p (the same in 
the two opposing lanes). The solid curves show the exact relation given by Eq. (8) whereas the broken 
curves above and below the curves for n = 3 and for n = 20 show the limiting relations given by Eqs. 
(8b) and (8a), respectively. Even for n = 3, one or the other of the limiting relations is a reasonable 
approximation over most of the range 0 < p < 1. 


. n(2 — p) | (1 — p) ] 
C, = we -. 
C2 (3 — 2p) , np(3 — 2p)(2 -- p)_|’ (8a) 





in which the second term in brackets is of order (np)~* « 1. The dominant factor is 
monotone increasing in p which means that the more left turns the better. This curious 
result comes about because in this range of p and n, the overpowering feature is the fact 
that the more cars there are in lane 2 that wish to turn left, the more chances there are 
for cars in lane 1 to turn left. In the limit p — 1, all cars turn left; there is no interference 
between the two lanes; and C, = C, = n. 

If n > 1 and p <1, particularly if np is of order 1 or less, the exponential terms in 
Eq. (8) are not small but by neglecting terms of order p and n~* as compared with 1 we 


find 


C,=C.~ n? + [1_— exp ( =i) (8b) 
3 Inp 


An interesting feature of this formula is the fact that C,/n, which, in a sense, represents 
the efficiency of the traffic flow, is a function of only np. In this range, C; is now a mono- 
tone decreasing function of np reflecting the fact that if left turns are rare, the oppor- 
tunities for left turns are also rare and the corresponding waiting period for a left turning 
car very long. For n > 1 > np (many cars per cycle but practically no left turns in 
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~] 


to 


one cycle), Eq. (8b) gives the proper limit C, ~ n but for n > np > 1 (small fraction 
of left turns but many left turns per cycle) both Eqs. (8a) and (8b) give C) ~ 2n/3. 

The most interesting feature of Eq. (8), however, is the fact that C,/n is a monotone 
decreasing function of n, a feature which would tend to favor small values of n. The 
competition that will arise between this effect and the obvious advantages of long traffic 
cycles suggests that under some conditions there might be an optimal cycle time. 

To investigate this possibility we note first that the traffic engineer has very little 
control over the values of p and p’ and so we consider these as fixed for any given traffic. 
One can, however, ch:nge n by changing the length of the green period of the light. 
This can be accomplished by changing the ratio of the green, yellow and red periods or 
by changing the time of a complete cycle. Changing the ratio of red to green times 
tends to favor one direction of traffic at the expense of the cross traffic. We consider 
here only the consequences of increasing the total cycle time 7 keeping the time of the 
yellow period fixed and the ratio of red to green periods fixed. Thus if ¢,, ¢, , and t, = yt, 
represent the times of the yellow, green and red periods respectively, 


T=t+0+y2, (9) 


we assume that ¢, and y are kept constant while ¢, and 7’ change. 

We may think of ¢, as a function of n; t,(n) is thus the green period that would allow 
n forward moving cars to pass through the intersection and 7'(n) is the corresponding 
cycle time under the above constraints. As a rough approximation, it is reasonable to 
assume that t,(n) is a linear function of n. For large n, after a steady flow has been reached, 
t,(n) will certainly be nearly proportional to n. The fact that a certain time is lost in 
reaching a steady flow can be roughly represented by a constant time loss. We therefore 


assume that 


t, = a+ bn. (10) 


The values of a and b depend upon the units of time, but a/b is dimensionless. We expect 
a/b to be a number of order 1. 
Substitution of Eq. (10) into Eq. (8) gives 
T(n) = (1+ y)b{n + [t, + (1 + ya)/b(01 + y)}. 
By measuring T' in units of (1 + )b, we can eliminate this factor so that 7 has the form 
T(n) =c+n, 

with 

e = [t, + (1 + yal/b(1 + ¥). (11) 


The constant c must be determined experimentally but we anticipate that c will be 
comparable with 1 and probably in the range of 1 to 3. 

The quantity of primary concern is the capacity of the intersection measured in 
cars per unit time, i.e. C,/7', and we wish to know if this has a maximum with respect 
to the parameter n. One can readily show from Eq. (8) that C,/T is a monotone increasing 


function of n for 


(Il — p) 
on =. — . ——— ( 1 
(2 — p)(3 — 2p)p 4) 


crm oO= 
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(the longer the cycle the better), whereas for c < c, , C,/T has a maximum for some 
finite value of n. There will be a maximum for sufficiently small p, and for c in the ex- 
pected range of 1 to 3 or so, the condition on p is 


XS 1/6c. 


We therefore anticipate that there will be a maximum C,/T if p is of order of about 
1/10 or less. It is not known to the author if this has ever been observed. 

The value of n at which this maximum occurs is the root of a transcendental equation 
and can best be found by plotting C,/7' or its derivative with respect to n vs. n. One 
can show, however, that the maximum occurs for np roughly of order 1. Figure 2 shows 
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FIG. 2 
Fic. 2. The existence of a maximum capacity C/T of cars per unit time as a function of n is illustrated 
by the curve c = 2 < co. The curve c = 4 > co is monotone increasing. The curves are shown as a func- 
tion of a continuous n although n is restricted to be integer valued. 


examples of a case in which C/T has a maximum and a case in which C/T is monotone 
increasing with n. 

Case II. p # p’, n > 1. Most of the qualitative conclusions deduced for p = p’ 
also apply for p + p’. The exact formula for N is quite complicated but for n > 1, it is 
possible to make approximations analogous to Eqs. (8a) and (8b). For n > 1 and either 
np > 1 or np’ > 1 or both, A} and dj are very small. The first term of Eq. (3) being 
multiplied by n dominates the other two terms and gives 


) 


1 

- Pp . n p’ — nN\/ 
vent 2a-p)+2a-phloa p’)/p\ 
ip’(1 — p)/p) 
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From this we obtain 


C, ee ©. 1! Sa |) (13) 


n° 1 + pil — 7), p +pi(l — p)/p 
and for C, a corresponding expression with p’ and p interchanged. These reduce to the 
leading term of Eq. (8a) if p = p’. 

We note that if p’ 0 (no left turns in lane 2) but np > 1 that C,/n—-0 
and C.,/n — 1. The significance of the latter is obvious; if there are no left turns in lane 2, 
there are no delays in lane 2. The fact that C,/n — 0 arises from the fact that with 
probability 1, only finitely many cars will pass the intersection in lane 1 before a left 
turning car appears and blocks this lane for the remainder of the cycle. Even for n > ~, 
we expect C, to remain finite and therefore C,/n — 0 forn > @. 

One can also see from Eq. (13) that C, is monotone increasing in p’ for all fixed 
values of p and monotone decreasing in p for all fixed values of p’. This is the type of 
behavior one expects since the more left turning cars in lane 2, the shorter the waiting 
time for a left turning car in lane 1 but the more the left turning cars in lane 1, the greater 
the number of delays. 

The other interesting limiting case for p ¥ p’ isn > 1, p< 1 and p’ < 1. By neglect- 
ing terms of order p’, p and n™"’, we obtain 














Ci __ + p’/p)_, {1 — exp [—(@ + p’ + pi’p*)nji p/p’) + (p'/p)""*] 
n° (1 + p'/p + p/P) 2n(p + p + pp”) + p'/p + p/P’) (14) 
4 {1 = exp [=@ + p’ = pp )\nI} (p/p) = (p'/p)"") 


nip +p —p p'’)\1l+p/pt p/p) 


, 


and for C,/n, the corresponding expression with p and p’ interchanged. For p = p’, 
this reduces to Eq. (8b). 

There are several limiting cases of this that are interesting. For n > np, np’ > 1 
(many left turns), 


oe St Ss 
n 1L+p/pt+p/p’ 


in agreement with the corresponding limit of Eq. (13), whereas for n > 1 > np, np’ 
(no left turns), C,/n ~ 1. For p/p’ < 1, we also find C,/n ~ 1 but for np of order 1 
and n > 1 > np’ (some left turns per cycle in lane 1 but none in lane 2) we obtain 


ig Re. (15) 
n np = 


This last relation obtains as a result of the fact that forn — © and p— 0 with np finite, 
there is still a certain fraction of the traffic cycle that will pass traffic before a left turning 
car blocks the lane. 

The question of existence of maximum capacities is much more complicated for 
p ~ p’ than for p = p’. Both C,/n and C./n appear to be monotone decreasing in n 
for all values of p and p’. There will certainly be a range of p and p’ values, with p and 
p’ sufficiently small, for which C,/T and C.,/T each has a maximum as a function of n. 
For p * p’, the values of n that give the maxima of C,/T' and C,/T will, however, be 
different. One would certainly not want to choose n less than the smaller or more than 





~I 


or 
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the larger of the two values of n defined by the two maxima, but any further specification 
of an optimal n would depend upon what function of the two capacities one wishes to 
maximize. 

Other models. The above model is neither the most difficult nor the simplest model 
of an intersection that one can investigate analytically. There are many relatively 
simple models for which one can evaluate at least the average flows when the queues 
are arbitrarily long. We consider below a few very simple examples. 

Suppose we have a multilane highway with at least two lanes of traffic moving in 
each direction, no left turn light and all left turn cars restricted to the left hand lanes. 
As long as all the queues remain non-empty, it is impossible to make a left turn at this 
type of intersection except at the end of the green cycle. The behaviors of the two oppos- 
ing traffics are independent of each other and can be studied separately. 

If we neglect right turning cars or assume that they leave the intersection at the 
same rate as forward moving cars, then in analogy with the previous model, we say that 
there shall be n cars leaving the intersection per cycle in each lane except in the one left 
turn lane. If a single left turning car will block the left hand lane, then j cars (j < n — 1) 
will pass the intersection per cycle in this lane if the first 7 — 1 cars in this lane are all 
moving forward but the jth car turns left after blocking the intersection for the remainder 
of the cycle; n cars will pass if the first m — 1 cars move forward, independent of what 
the nth car does. 

If left turning cars occur in the left hand lane with probability p, then the prob- 
ability of j cars passing the intersection in this lane is 


p1—p)’' for j<n-—1, (l1—p)"" for j=n. 
The average number of cars to leave in this lane per cycle is 
p >, jl — p)i* + n(1 — p)”" = pf — (1 — p)’). 
If there are m lanes of traffic moving in the same direction, the total capacity of the 
intersection in this direction is 
C=(m—1)n+p [1 — (1 — p)"). (16) 
If, instead of having just one left turn car block the intersection, we assume that 
there is a space for one left turn car to wait but that a second left turn car will block 
the left hand lane, then the probability of 7 cars passing the intersection per cycle in 
this lane is 
(j — 1p? — p)*”? for 2<j n-1 
(1—p)"'+(m—1D0 — p)"*? for j=n. 


The average number of cars to leave in this lane per cycle is 
P > jG — DA — p)*? +n — pp)” + n(n — Ip — p)"” 

= 2p"[1 — (1 — p)"] — n(l — p)*”. 
The total capacity analogous to Eq. (16) is 


C = (m — 1)n + 2p"[1 — (1 — p)"}] — n(l — p)*”. (17) 
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Each of the above situations can lead to an optimal cycle time in the sense described 
previously for the two lane highway. If we take 7'(n) as in Eq. (11), then C/T from 
Eq. (16) has a maximum as a function of n if cp(m — 1) < 1 and C/T from Eq. (17) 
has a maximum if cp(m — 1) < 2. 

One could also extend these calculations for the multilane highway to include some 
effects of right turning cars. It would be reasonable to assume that right turn cars appear 
only in the right hand lanes and that the capacity of this lane is independent of what 
happens in other lanes. This problem can therefore be considered separately from the 


above. 
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THE THICKNESS OF CYLINDRICAL SHOCKS AND THE PLK METHOD* 


BY 
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Summary. Cylindrical shocks occur when a viscous heat-conducting gas flows 
radially in a plane. This is a singular perturbation problem in which the perturbing 
parameter is the reciprocal of R, , the Reynolds number of the flow. It is shown that for 
both inward facing shocks (source flow) and outward facing shocks (sink flow) the 
shock thickness is of order R7'S~* log (R,S*) where S is the shock strength. This is 
contrary to results for sink flow which have been obtained by the use of Lighthill’s 
technique for rendering approximations uniformly valid—the PLK method. It is shown 
that this method fails when applied to singular perturbation problems of the type 
discussed here in which the small parameter multiplies the highest derivatives. 

Introduction. The steady radial two dimensional flow of a viscous heat conducting 
gas has attracted some attention in recent years [1, 2, 3]. Perhaps the main interest 
has followed from the fact that if an inviscid gas with zero thermal conductivity existed, 
it could only flow in such a manner exterior to the sonic circle which is a limiting line 
of the flow. Outside the sonic circle the speed could either be supersonic everywhere 
to infinity where the vacuum speed is attained or subsonic everywhere to infinity 
where the speed is zero. It is apparent that this is a singular perturbation problem 
because the flow of a real gas in the limit of vanishing viscosity and thermal conductivity 
should not have this character. 

Since there is only one independent variable for this simple configuration a treatment 
of the problem is feasible. A Reynolds number R, may be defined in terms of the mass 
flow and the viscosity, and the solutions sought when R, is large and the Prandtl number, 
a, is finite. 

It has been shown by the writer [2] that when the gas flows outwards, source flow, 
the limit of the solutions as R, ~ © exhibit no ambiguities. Although there is a singular 
stagnation circle, the ‘source’, which falls outside the scope of the Navier-Stokes 
equations, the important feature is that a typical solution curve follows the supersonic 
branch of the inviscid solution outwards from the sonic circle for some distance (dependent 
on external boundary conditions), then crosses via an inward facing cylindrical shock 
to the subsonic branch of the inviscid solution to complete its journey to infinity. If R, 
is large but finite the shock is not a discontinuity and has a thickness of order 
R-'S~* log (R,S*), where S is the shock strength. 

More recently Wu [3] has investigated the case of inward flow, sink flow, and finds 
an analogous situation. In the limit R, — ©, a typical solution curve follows the super- 
sonic branch of the inviscid solution from infinity in for some distance, then crosses via 
an outward facing cylindrical shock to the subsonic branch of the inviscid solution. 
The continuation beyond the sonic circle exhibits a singular vacuum circle, the ‘‘sink.” 
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However he finds that for finite R, an outward facing shock can only be of strength of 
order R;*”* and have a thickness of order R,*”’ 

Although it would be expected that the results should agree, a proof seems to demand 
a discussion of the whole flow field. It is for this reason that the difference is important, 
for instead of the topological methods used by the writer, Wu has employed an analytical 
method, the PLK method, to investigate the flow field, and this approach has been 
endorsed by Tsien [4]. Hence the procedure which will be adopted is to show that these 
results are wrong in general, Sec. 2, and to indicate the reasons why invalid results 
were obtained, Sec. 3. 

The governing equations for the flow are found in suitable form in Sec. 1 and after 
a brief consideration of the singularities for the general case, the discussion centres on 
the iso-energetic flows, Sec. 2, for a particular value of o close to ?. For this case a topo- 
logical discussion of the governing equations is possible, and as has been pointed out 
by the writer [2] variations in o will not affect the nature of the results. After a brief 
review of source flow, entirely analogous results are derived for sink flow. In particular 
it is shown that, generally, the outward facing cylindrical shocks which occur have 
a thickness, A, of order R;*\S~* log (R, S*), while, in their interior, the maximum non- 
dimensional velocity gradient, V, is of order R, S*. It is only for the very weak shocks 
for which S is of order R7’” that Wu’s results, A = O(R7**), V = O(R7”), hold, 
and these are still covered by the general result. An interesting sidelight here is that 
the weakest cylindrical shocks which can occur are of strength of order R7'”’. 

Now Wu’s results are obtained from similarity solutions valid for speeds near sonic, 
but the essential feature is that the PLK method is used to show that all solution curves 
join on to these similarity solutions. Hence it is his application of this method, the 
Lighthill ‘‘technique for rendering approximations uniformly valid’”’ [5], which comes 
under fire in Sec. 3. 

In Lighthill’s expository paper this technique was applied to singular perturbation 
problems which arise because singularities are wrongly placed by the zeroth perturbation. 
There is no reason a priori to expect that it will be of any use for the different class of 
perturbation problems which are singular essentially because the highest order deriva- 
tives are dropped in the equations governing the zeroth perturbation, as in our case. 
The effect of the highest derivative is felt almost everywhere in this problem, and neither 
a conventional perturbation expansion nor a Lighthill perturbation expansion, which 
still loses the highest derivative in each perturbation, leads to valid solutions. This is 
brought out clearly in Sec. 3. A critical examination of the application of the PLK method 
to the sink flow problem is followed by an example of its application to a simpler equation 
of the same broad type, for which the explicit solution is known. The manner in which 
it fails to provide any improvement to a conventional perturbation expansion is then 
easily seen. 

The moral, of course, is that it is useless to expect the PLK method to be a panacea 
for all singular perturbation troubles. Their nature must be assessed as suitable for the 
treatment first, and it is possible that some other applications of the method should 
be critically reexamined. 

1. The fundamental equations. 


lor purely radial two dimensional flow of a fluid 
the speed, u, is only a function of r, the radial distance from the origin. Let p, p, T, u, 


o, R, c, and a denote respectively the pressure, density, temperature, viscosity, Prandtl 
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number, gas constant, specific heat at constant pressure and local speed of sound. Then 
the equations for the conservation of mass, momentum and energy take the form [2] 


pur =x, (a constant), (1.1) 
du, dp_ 4 a ( du) 2udu _ Quu _ 2d (“*) 
pn ar + dr 3adr\" adr + r dr r° 3dr\r/’ (1.2) 
and 
lL ss 2 du du c,ur dT 
ry ba 2 = — <= @ ee i ee ee i € 
c(e,1 4 5 oo 3 sul Tp a u) 2uur 7 Age kI, . (1.3) 


In addition, if the fluid is an ideal gas the thermodynamic variables are related by the 
equation of state 

p = RoT. (1.4) 
By analogy with inviscid flow the constant J, is essentially positive while «x is positive 


(negative) for source (sink) flow. 
In terms of the non-dimensional variables 











w= |u| (2e/,)7"'”, 6=c,T/I,, a = 4u(1 + 8)/(3x), (1.5) 
with 
8 = (vy — Il)/v + D = R/(2c, — R), (1.6) 
elimination of p and p with the aid of (1.1) and (1.4) yields the two equations 
1+6dw,d (£) _ (f Ldw _ ¥) (#2 m 2) da 
r dr + dr\wr) — “\ dr + rdr r sel dr 2r/ dr’ (1.7) 
and 
( a - ___3ar dé = 2aB8wrdw _ 
iad. |, + Bay “)u 411+ Bo dr 1+8dr — i (1.8) 


We may regard | a@ | as the reciprocal of the Reynolds number of the flow. 
If a is placed equal to zero, these equations together yield the inviscid approximation 


—r(dw/dr) = w(1 — Bw*)(1 — w’)™ (1.9) 
in which w = 1 is the sonic speed and w = 6°‘ is the vacuum speed. The solution of 
this equation is 

(r/r.) = w [1 — Bw?)/ — ayo“ (1.10) 


where r, is the radius of the sonic circle, and shows that external to the sonic circle the 
field is doubly covered, so that at each point the speed may be either subsonic or super- 
sonic, while inside the sonic circle there is no solution. When r — © either w — 0 or 
w — 8 '” (see Fig. 4). Hereafter we shall suppose r to be rendered non-dimensional in 
such a manner that the radius of the sonic circle is unity in this corresponding inviscid 
flow—this does not alter the form of the equations (1.7) and (1.8). 

l’or simplicity, we shall only consider the flows for which a is constant. This implies 
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in particular that the viscosity is independent of the temperature but previous results [2] 
suggest that this will not affect the general nature of the conclusions to be drawn. In 
this case, if o takes the particular value o, where 


40, = 3(1 + a/(8 + 1)], (1.11) 
it is possible to obtain an energy integral. For if we define 
E=6-—1+ All + a/(6 + 1)Jw’, (1.12) 
so that 1 + Z£ is essentially the total energy of the flow, Eq. (1.8) yields 
E = cpt (i tbiee/ (3a) (1 13) 


where C is an arbitrary constant. 
Now, consider the equations obtained in terms of the new variable 


V = -r dw/dr, (1.14) 
so that V is essentially a velocity gradient. They are 
aw’ V(dV/dw) = aw* — (1 + B)w*°V — wV(d0/dw) + Vo — we (1.15) 


and 


(1 + B)(6@ — 1) + BU + B + @)w* + 8a(4c)"'V(d0/dw) + 2o6wV = 0. (1.16) 
It may be shown that the singular points of these equations are at w* = w? = [8 + 
a(1 + 28)/(1 + 8)]"’ where 6 = 0. = aw? and V = 0, and at w = 0 where @ = 1 and 


V = 0. 
In the w-V plane the point w = 0, V = 0 isa double point*, and the possible gradients 
= 1. For both a > 0 (source flow) and a < 0 


there are (dV/dw) = ~ and (dV/dw) = 
(sink flow), there is only one solution curve of infinite slope through the origin, namely 


w = 0. For source flow there is an infinite number of solution curves with (dV/dw) = 1 
through the origin, but there is only one solution curve through the origin with this 
slope for a < 0. 

The singular point at w = w, = 8°” (the vacuum speed), V = 0, is a triple point 
in general in the sense used above. For the case of source flow, there is one curve through 
this point with slope (dV/dw) = X, , on which (d6/dw) = yu, , corresponding to r — © 
(this is the physically interesting one and corresponds to one branch of the inviscid 
solution) while corresponding to r — 0 there are curves described by 


V~ rw — w) + C | w —w, Pe, | (1.17) 
6— 0. ~ us(w — w.) + C[(u2 — us)/(2 — As)] | w — wv pene) 
and if C ~ 0, 





V~ Aa(w — We); = 6. —_ b2(w ies W.); (1.18) 
where 


A; 
ho = —(1+ Bla’ +O0(1), me = 16c8a[3(1 + 8) — 40(1 — A) + Oe’) 


3A; = —4o(1 — Bla’ + O(1), 3us = [40(1 — 8) — 3(1 + B)]w. + Of). 


(28)/1 — 8) + O@), mw = —26w. + Of), 


yp (1.19) 
7 


*This is used in the sense that there are two possible solution curve gradients through V = w = 0. 
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If ¢ = oo, then it may be shown that on these curves 


E sia ional 


but from (1.13), in this case 
E tis cr thle 

for all r, which is unacceptable if the inviscid solution is to be a limiting case when 
a — 0 for large r. Hence C = 0, Z = 0, and the only curves left are 

V~ 3(w = W.), (r — 0), 

V ~ A, (w — w,), (r— ~). 
On the other hand, for sink flow, this point is again a triple point, but all the curves 
through it correspond tor — o, If 


40 < 3(1 + B)/(1 — 8) 


they are 





V ~ ,(w — w,) + (As — AK, | w — w, | + A, — ADK, | w — w, |, 
with (1.20) 
0— 0. ~ u,(w — w.) + (us — mk, | w — w, [+ (ue — maka | w — es 
(which corresponds to the supersonic branch of the inviscid curve), 
V ~As(w — w.) + (Ae — As)Ch, | w — w, |”, 
with (1.21) 
6 — 0, ~ us(w — w.) + (ue — ws)Ck, | w — w, |, 
and, if k, ¥ 0, 
V ~ d.(w — w,) 
(1.22) 


with 
6— 0. ~ u(w — w,). 
When ¢ = oo , Hs = m , and we find that 
E w~ Ck Ot? "'*' for ro. 


Again from (1.13) we have 


E _— Cr” Otel 
4 


for all r for which the solution exists and it is tempting to argue that in order to keep E 
bounded when r is small and | a | — 0 we should put C = 0 and hence k, = 0. But in 
the absence of a general solution it is by no means certain that r is ever less than one on a 
solution curve yielding a shock, in fact in the isoenergetic case it is not true anywhere 
where the solution has physical meaning. 

Thus, while for source flow the special case ¢ = a» yields the result that all flows of 
bounded total energy have constant total energy it is possible that for sink flow all 
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flows with ¢ = o» have bounded total energy. In order to achieve any simplification 
then we will consider only the sink flows with constant total energy and it seems plausible 
that the results obtained for this case will be indicative, at least, of those holding for 
neighbouring flows—and anyway, they will be sufficient for one of the main purposes 
of this paper. 

For this particular case, the singular point at w = w, becomes a double point, through 
which the solution curves are 


V ~r,(w — w.) + As — ADK, | w — w, | (1.23) 
and 
V ~ i,(w — w,) (1.24) 
on both of which 
6— 06.~ u,(w — w,) (1.25) 
since now k, = 0 and wz = uw, when o = op. 


Note that if the isoenergetic assumption is not made, then from (1.20) there is a 
double infinity of curves (involving two arbitrary constants k, and k,) with slope i, 
which therefore are asymptotic to the inviscid curve. 

2. The isoenergetic flows with o = o,. We have now the relation 


6=1— sl +a/(8B + 1)]w’ (3%) 
and this may be used to eliminate 6 from (1.15) to yield the first order equation 
aw V(dV/dw) = V{1 — [1 — a8/(B + 1)Jw’} (2.9) 
— wil — [8 + a(l + 26)/(1 + B)Ju’}. 


There is only a trivial loss of accuracy in replacing the coefficients of w* in the brackets 
by 1 and 8 respectively [2], and the equation which will be considered subsequently 


becomes 
aw’ V(dV/dw) = V(1 — w’) — w(l — Bw’). (2.3) 
If a is placed equal to zero here we obtain the inviscid equation 


V = w(l — Bw’)(1 — w’)"', (2.4) 


another form of (1,9). 

The equation (2.3), although not soluble explicitly, lends itself to a topological 
discussion of the solution curves which has been carried out by the writer [2] for the 
case of source flow. The solution curves of the equation run as in Fig. 1 and exhibit the 
following features. Those of physical interest ‘start’ from the point w = V = 0 (a stag- 
nation point at r = «). A typical curve lies close to the subsonic branch of the inviscid 
curve up to some subsonic speed w, say, then rises steeply to a large positive value of V 
and returns to the neighbourhood of the supersonic branch of the inviscid curve near 
w = w,' which it then follows nearly back to w = 1. Ultimately the solution curve 
becomes asymptotic to w = 0 with V — — ~. In the physical plane the steeply humped 
portion of the curve corresponds to a shock of small, but finite, thickness while the part 
for which w — 0 with V — — © corresponds to a singularity lying close to r 1 from 
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Fig. 1. Source flow: Solution curves in the w-V plane 


which the fluid issues with infinite density. Thus in the region of interest we have effec- 
tively parts of the inviscid solution branches joined by a shock. 
It was shown [2] that for such a typical solution curve, if 


(1 — w,)™" = ofa’) (2.5) 
that is, the shock strength, 1 — w, , is of larger order than a'* then 
A = Ofa(1 — w,)™ log [(1 — w,)*a™']] (2.6) 
where A is the shock thickness. Furthermore, in the interior of the shock 
Vena = O[(1 — w,)’a™’). (2.7) 
The case 
1 — w, = Oe”) (2.8) 
was not discussed, but a similarity argument may then be applied and yields the results 
A = O(a”), (2.9) 
Vinx = O(0"”), (2.10) 


consistent with (2.6) and (2.7). 
lor sink flow, the governing equation is written 


la | wV(dV/dw) = —V(I — w’) + wil — Bu") (2.11) 
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and this has been discussed in some detail by Wu [3]. The curve of zero slope is the 
inviscid solution curve, €, , given by 
V = w(l — Bu’)(1 — w’)" (2.12) 
and the curve of inflexions, @, , has the equation 
QV? — wil + Bw) V? + Jal’ (1 — Bw®)[V1 — w’) — w(l — Bw’)] = 0. (2.13) 


Provided that 
11 —w| = oa’) (2.14) 
€, has branches given by 
V = w(l — Bw’)(1 — w’)" fl + o(1)] (2.15) 
and 
V = 4(1 — Bu’)'?(w? — 1)'7(2 | @ |)" [1 + 0(1)], w>l; (2.16) 
while it has infinite slope at 
wr~rl1+3.2°“4(1 — Bp)? la |”, V~ -(1 — p74] a)”, (2.17) 
and crosses the line w = 1 at 


V 











Fig. 2. Sink flow: Solution curves in the w-V plane 
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V~(l — p72 /al)~™”. (2.18) 


Hence the curve runs as shown in Fig. 2. It is easily verified from power series expan- 
sions that although the infinite family of solution curves passing through w = 67” 
have the same slope 28/(1 — 8) [see (1.20)] as @, and the branch (2.15) of @, , they 
lie initially below both those curves. Thus the solution curves may be sketched into 
the figure. 

Now consider a typical solution curve from the point w = 6°”, V = 0, which 
passes through the point P(b, c) on the lowest supersonic branch of the curve of inflexions, 
as in Fig. 3, where 


[b-1|* =olje|™), (2.19) 











| 
| 
| 
| 
l 
| 
A A 





Fig. 3. Sink flow: A typical shock curve in the w-V plane 
so that from (2.16) 
c= —(1 — pb*)'°(b? — 1)'°(2 | a |) *7[1 + o(1)). (2.20) 


Since P(}, c) is a point of inflexion then the solution curve certainly lies to the right of 
its tangent there, 


V=c+(b’ — 1b? lal" [1+ 011 — Bh%)c'(h — 1)"Jw—d, (2.21) 
for w > b, and to the lefv of this tangent for w < 6. On the other hand, on the line PQ, 
V=ct+(b — 1b* lal" [1 + b(1 — Bb*)c"(b* — 1)" ]1 — dw — By), 

S<e< f, 


(2.22) 
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it may be shown that the slope of the solution curves is greater than the slope of the 
line for 

Wo = b— 3 'Bbe(b° — 1) Cw b, (2.23) 


hence the solution curve through P cannot cross it in this range and must lie to the right 
of PQ. 
A repeated application of this kind of reasoning shows that (see lig. 3) 
(i) the curve lies to the left of the line PR, where 


MR = 2°*0°7(1 — Bb’)**(b? — 1)'"9(b)"” | a |" (2.24) 

and 
(RL/ST) = 2°°0'7(1 — Bb?)-'*(B? — 1)7*“9(b)'? | a |", (2.25) 

where 

g(b) = 1+ (1 — 36)b’ + Bb‘ > O, (2.26) 

and 
ST = b(1 — gb’)(b? — 1)"; (2.27) 

(ii) the solution curve crosses the w axis at w = d, where 

1—d=0(b—-)), (2.28) 


and lies to the left of the curve P’Q’ which has the equation 


V = -2(1 — pd)'7d-"? | a |"? (w — od)” (2.29) 
— 23°'(1 — d’)d* |a|"' (w — d), 

with 

We —~d=b'd(l—d); (2.30) 
(iii) the solution curve lies to the left of the line P’R’ where 

SR’ = d(1 — Bd’)(1 — d’) "(ad)"? | a |'”, (2.31) 
(L’R’/P’S’) = g(d)'7 01 — d’)*? la |”, (2.32) 

and 
P’S’ = d(1 — Bd’)(1 — a’). (2.33) 


Hence the solution curve through P lies close to the inviscid curve for w slightly 
greater than b, from (2.24) and (2.25), and for w slightly greater than d(V > 0), from 
(2.31) and (2.32), and between these points, on the whole, the velocity gradient is large. 


In fact, the curve in this range represents a shock, and the solution in the r-w plane 


is shown in Fig. 4. 
We may estimate the shock thickness by means of the bounding curves. For con- 
venience, we define the shock thickness as the distance between P and P’. Since 


V = —rdw/dr 
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Fig. 4. Sink flow: A typical solution curve in the r-w plane 





the preceding results show that log (rp/rp-) is bounded above by 
Q’ . P 
/ tan dw, +* [ | V | dw + (d — b)[min. (| Vo. |, |] Ve I", 
Pp’ iJo 


and below by 
Ps 
(1 —d)|V, |’ + | | V |7* dw. 
A 


It is easily verified that the dominant term in each of these bounds, and hence log (r,/rp-)*, 
is of order 
(b — 1)7' | a | log [(b — 1)°(1 — Bb)" | al7'] (2.34) 
since (1 — d) is at least of order (b — 1). 
Now frp will differ little from the value of r at w = b on the supersonic branch of the 
inviscid curve, namely 


b-"[(1 — pb’)/ — ayo", (2.35) 


so that finally the shock thickness 


‘ 1 | ' Q},2\ — (1-B)/ (28) 
A fp: — Tp = O}(b — |) feid-— Bb’) (2.36) 


‘log [(b — 1)? | a |~* (1 — Bb?)"'}}. 


This result may be used to show that Eq. (1.7) integrated between rp and rp, yields 


2 


bd~1+ fal’? 2°70 — gb’*)'7(b? — 10°” (2.37) 


which is the Prandtl relation to order | a |'” 


From (2.21), (2.23) and (2.30) the maximum value of | V_| attained within the shock 
lies between | V, | and the minimum of | Vg | and | VQ. |, which are all of order 
(b— 1)? lal. (2.38) 


*To extend this result to the neighbourhood of the inviscid curve near the points Z and L’, more 
precise bounds than those given in (i) and (iii) are needed, but the result (2.34) still stands. 
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But we note that at P, 


1 V| =Ol(b — 1)? |a|""7] (2.39) 





from (2.20), and hence the shock curve cannot be obtained by a similarity argument. 


However, if 


(b— 1) = O([a|"”’) (2.40) 


the above results suggest that 


|V| =Ola|”’) (2.41) 


throughout, and indeed, with the change of variables 


V=l|ea ace v, (2.42) 
' ! 3 ‘ € 
w-l=lal'’*a, (2.43) 


similarity solutions of (2.3) may be obtained [3] to yield shock type curves for which 
A = O(|a |”). (2.44) 


But this is obviously a very restricted case, and it is important to investigate the 
reason why these were the only solutions obtained by the PLK method. 

3. The PLK method. When the governing equations (including boundary con- 
ditions) of a problem contain a small parameter, say a, the standard perturbation 
technique is to expand each dependent variable as a series in a. The equations for the 
vth perturbation are obtained from the coefficients of a’ in the governing equations— 
in particular, the zeroth perturbation satisfies the unperturbed equations. 

These perturbations may exhibit, at a point (or on a curve), a singularity which 
becomes worse as the order », increases. If this occurs when there is only a ‘slight’ differ- 
ence between the perturbed and unperturbed equations and in particular, they are of 
the same order, Lighthill [5] has shown that in many cases a uniformly valid solution* 
may be obtained as follows. Expand also appropriate independent variables as series 
in a in terms of new independent variables. The coefficients of a’, except for a’, are 
unknown functions and are used to cancel the worst effect of the singularity, so that in 
terms of the new independent variables each perturbation is no more singular than the 
preceding one. A transformation back to the original independent variables yields the 
uniformly valid solution. 

The singular perturbation problems just described are of this nature essentially 
because a singularity is incorrectly placed by the unperturbed governing equations 
(it may even lie on another sheet of a Riemann surface), and the Lighthill technique 
(the PLK method) places the singularity in its correct position. The problem which 
concerns us is of the singular perturbation type because the unperturbed equations are 
of lower order than the perturbed equations, and are singular at w = 1 whereas the 
perturbed equations are singular only at w = 0 and near w = 6’”’, see, for example, 
Eqs. (2.3) and (2.4). 

It is this spurious singularity at w = 1 which gives rise to a superficial resemblance 


*Though not necessarily to all orders, see Tsien [4], p. 334. 
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to the problems treated by Lighthill. An attempt to solve Eqs. (1.7) and (1.8)* by means 
of conventional perturbation expansions shows that successive perturbations are more 
and more singular at w = 1. Hence, a new variable [3], ¢, is introduced and the expansions 


wt) =§E+]al we) + lal? we) +---, (3.1) 
logr = not) + | a| m) + [a |? me +--+: , (3.2) 
6(E) a 6o(é) + la | 6,(&) + la |? 6.(E) + cre (3.3) 


are assumed. 

Upon substitution into the differential equations (1.7) and (1.8) and equating co- 
efficients of powers of | a |, two differential equations which must be satisfied by w, , 
n, and 6, are obtained, for each r. In particular the zero order equations are 

(1 + BE + £05 — A — ENO = O, (3.4) 
6.—1+ p88 =0, (3.5) 
provided that »{(£) # 0. The solution for 7 is 
no(é) = —log & — (28)"'(1 — B) log [(1 — A€)/(1 — 8)], (3.6) 
where the constant of integration is chosen so that no(w) is the inviscid solution for log r. 
The »,(£) are now determined in turn so that the w,() have the lowest possible order of 
singularity at ¢ = 1**. It is found that near ¢ = 1, 

mé) =O(1-"), mm) =Ol1-H“*)], nm@ =O -—H***). (3.7) 

Since it is noted that the singularity near w = 8~'”” cannot be accounted for by the 
expansion procedure adopted, the series (3.2) is finally written 

logr = nt) + | a] m&) +++ $C] 1 — pe [OP's (3.8) 
where the last term is obtained from a discussion of the singularities of the isoenergetic 


equation. 


° ° ° ° 1 . 
It is suggested that this expansion is convergent when | 1 — ~! >| a@ |’, since 
successive terms are of smaller order in | a |, and moreover, that it remains convergent 
even when | 1 — —! = k| @ |, where k is indeed of order unity, although only the first 


few terms are known explicitly. 

Now, it may be verified, at least for the isoenergetic case, that the series obtained 
are divergent, even when | 1 — ~| >| a |'”*, but the important points can be brought 
out much more clearly by the examination of a simpler equation which, while it exhibits 
the main features of (2.11), possesses an explicit solution. Consider the equation 

a(dy/dxr) = —zy — 1 (3.9) 
in which a > 0 is a small parameter. The equation has no singularities in the finite 


plane, but the equation governing the zeroth perturbation, 


zy +1=0, (3.10) 


*The formal development which follows is essentially the same for the general case and the iso- 
energetic case, except for the term with arbitrary constant in (3.8). 

**It would seem equally plausible to determine the w, so that the », are no more singular than 
no but no essential difference results. 











90 H. C. LEVEY [Vol. XVII, No. 1 


obtained by placing a equal to zero, is an equation of lower order than (3.9) and is 














singular at x = 0. 

The solution curves of (3.9) are sketched in Fig. 5 and it is seen that for y < 0 thei 
general character resembles those of (2.11). The point z = © corresponds to w = 6 '”* 
and the curve 

y = —" (3.11) 

corresponds to the inviscid curve (2.12). 

x 
Fia, 5. The solution curves of ay’ = —zy — 1, 
An attempt to solve (3.9) formally by means of an expansion of the form 
y = yo(x) + ay,(x) + --- (3.12) 
leads to 

y= —2 — ag? — ees — 1-3-5 +++ (Qn — laa" «ee, (3.13) 


which is divergent, successive terms being more and more singular at x = 0. The exact 


solution of the equation for the curve which passes through the point 2% , Y is 


y = Yo exp [(r5 — 2”)/(2a)] — a” exp [—2"/(2a)] | exp [t’/(2a)] dt, (3.14) 
and it may be verified that (3.13) is the asymptotie expansion of (3.14) for x large enough, 
yo , but it is not a valid expansion for any solution for all x > 0. In particular 
1— wt! = O(! a |’) before), all of the terms in (3.13) 
'/3) but, it can be verified 


for all x5 , 

1/2 
when z = O(a ’*) (note that 
are O(a ’’*) [note that all terms in (3.2) are of order | a 
from (3.14), the remainder term is also O(a ‘’*) and becomes infinite with r for any 


Zo , Yo . Thus the series (3.13) is never an adequate approximation to a solution in the 
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critical region when x = O(a'’’). (The addition of a term [exp (—2*/(2a)] does not 


improve matters). 
Now apply the PLK method to this equation. A new independent variable is chosen 


so that 
x= + az,(t) + a’x,(t) + --- (3.15) 
and it is assumed that 
y = yot) + ays) + a*y.(f) +--+. (3.16) 
The x,(£) are to be determined successively so that the y,(€) are no more singular than 


yo(£) at € = 0. The substitution of these expansions into (3.9) yields the equations 


ty +1 = 0, (3.17) 
(EYo + I)a{ + 21Yo + fy, = — Yo ; (3.18) 
(Eyo + I)az + (iyo + Em )zi + EY2 + MY: + L2Yo = —Yi,***- (3.19) 


[It is possible to annihilate all of the y, for n > 1 by suitable choice of the x, to obtain 
the formal solution 


y= —¢", (3.20) 
em ttof’ + ab +--+ + anger’ +---, (3.21) 
where 


dom _ (4m pe) 2 A2m—1—0%e41 + to: |, m > 2, 


— (3.22) 
m-—1 
Gens, = 4m > a ee m> 1. 
It follows from the last relations that 
a, > 2” *(n — 1)!, n> 2, (3.23) 


so that the series (3.21) is certainly divergent, in fact, ‘‘more’’ divergent than (3.13), 
and the successive terms are more and more singular at £ = 0. The series does not enable 
the neighbourhood « = O(a'”*) to be approached, and all that has been achieved by 
the PLK method is the exchange of an invalid expansion of one variable for an invalid 
expansion of the other. The crucial point is that the small parameter multiplies the 
highest derivatives in the differential equation. Therefore the highest derivatives are 
dropped in the sequences of equations set up both by the conventional perturbation 
analysis and by the PLK method; this is illustrated in (3.17) to (3.19). But the explicit 
solution (3.14) shows that the highest derivative in (3.9), despite its small coefficient, 
is at least as important as any other term over a large part of the field. 

For the cylindrical shock, the situation is similar. In Wu’s application [3] of the PLK 
method, the highest derivatives are dropped from the sequences of equations, but their 
importance is indicated by the topological results, at least for the isoenergetic case, and 
it is not surprising, therefore, that the expansions diverge before a neighbourhood of 

| is reached. 
The term in (3.8) with the arbitrary constant serves to represent the effect of the 


c 
S 








92 H. C. LEVEY [Vol. XVII, No. 1 


highest derivatives in a neighbourhood of — = 8” '”* for the isoenergetic case [the discus- 

sion of the real singularity near w = 6" ‘”’ in Sec. 1 shows that near ¢ = 8~'””, in general, 
Y | 2 |As/d ait hake ° ° 

log r ~ m(é) + Cll — BE I" + D1 — BF | ], but even with this term, the 


expansions are only asymptotic to the solution in that neighbourhood, since the series 
diverge just the same. 

Seeing that all the terms 7,(&) in the (invalid) expansion (3.8) are 0(| a 
'1 —¢| = 0(! a |’), Wu [3] concludes that for all cylindrical shocks, log r = O(! a |?”*) 


> 


3 
) when 





when ] 1 — €| = 0(| a |*”*), and he constructs similarity solutions 
logr = | a |** &, (3.24) 
w= 1+ [al wie) +---, (3.25) 
@=1+ fal’? oxe*) +---, (3.26) 


the validity of which is not contested, in the range 


a , 





l1—kla|[’<t<1+k 
The similarity solutions exhibit a shock-like character and imply that across such shocks 
Aw =Ala|'” 


and 
Alogr=Bla|™”. 


To determine typical values of the constants A and B, Wu [3] takes the divergent 
expansion (3.2) to represent the particular solution corresponding to C = 0 in (3.8), 
and uses the few terms, which are available, of it and of the corresponding expansions 


for w and @ to compute values of log r, wand @at& = 1+ k| a |'” to serve as boundary 


values of the similarity solution (3.24) to (3.26). 

It is understandable that such a procedure can only lead to the selection of a restricted 
solution class defined by the similarity assumption. What follows from the existence of 
the transonic similarity solutions is that shocks of strength of order | a |'”* have a thick- 
*/° This restricted class corresponds to the solution curves in Fig. 2 


ness of order | a | 
'? of its tip. 


which pass through the curve of inflexions within a distance of order | a | 
For the stronger shocks the results (2.36) and (2.38) hold. 

It is interesting to note that it is impossible for a cylindrical shock of zero strength 
to exist—at least for the isoenergetic case—for there is a continuous transition from 
weak shocks with strength 0(| a |'’*) to ‘incomplete’ supersonic (subsonic) ‘shocks’ of 
strength O(| a |'”*) for sink (source) flow as the point P, Fig. 3, moves along the curve 
of inflexions to the tip and back again on the other branch. Thus the weakest shocks 
obtained in cylindrical flow have a strength of order | a@ |'”° 

4. Conclusion. It has been shown, in agreement with the result for inward-facing 
cylindrical shocks 
at least for the isoenergetic case, have a thickness given by 


(source flow), that outward facing cylindrical shocks, (sink flow), 


A = O[S* | a| log (S* | a |~’)], (4.1) 
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where S is the shock strength and | a | may be regarded as the reciprocal of the Reynolds 
number of the flow. The maximum nondimensional velocity gradient within the shock is 


(dw/dr) max = O(S’? | a |~'). (4.2) 
The weakest cylindrical shocks which can occur are of strength 
S = O(| a |'”) (4.3) 
and then 
A= O(|a|’"), (4.4) 
(dw/dr) max = O(| @ |”), (4.5) 


in agreement with the results for arbitrary strength. It is only in this case that a similarity 
solution can describe the shock interior. 

Singular perturbation problems of this type cannot be treated by a straightforward 
application of the PLK method. This approach,leads to the invalid result that all cylin- 
drical shocks have a strength of order | a |'”*. 

Acknowledgement. Acknowledgement is made to the Chief Scientist, Department 
of Supply, Australia, for permission to publish this paper, and to Professor R. E. Meyer 
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BOOK REVIEWS 


(Continued from p. 66) 


placed at the end of the volume. Chapter I is entitled “Volterra equations”. Topics treated include 
Volterra’s integral equation of the second kind: f(r) — Jt K(z, y)f(y) dy = g(x); Volterra’s integral 
equation of the first kind: fi K(a, y)f(y) dy = g(x); relationship between Volterra integral equations 
and linear differential equations of the type 

d"u ‘8 


1 

u . 
dx” az t+ es + .4,(z)u = F(z); 
ax 


a 
+ a,(x) 
integral equations of the Faltung (‘‘closed cycle’’) type, where the kernel K(z, y) = k(x 


dx 
- y); Volterra 
= F(x, y)/(z — y)*, and non-linear Volterra equations. Of 


equations with kernels of the type K(z, y) 
particular importance are the applications to the asymptotic evaluation of Bessel functions for large 
values of the argument and also to the more difficult asymptotic evaluation of Bessel functions whose 
order and whose argument are almost equal. Chapter II bears the title: “Fredholm equations.” It 
begins by the solution of the Fredholm integral equation of the second kind; f(z) — A fo K(z, y)f(y) dy = 
g(x), which leads to Neumann’s series; together with the application to a particularly simple example. 
The Fredholm theorem for general ZL; kernels is proved by using E. Schmidt’s idea of using the Neumann 
series to pass from “degenerate” (Pincherle-Goursat) kernels K(x, y) = )-%e1 Xi(z)¥x(y) to the general 
kernel. Next follow Fredholm’s formulas and a section on the numerical solution of integral equations. 
The chapter concludes with the Fredholm solution of the Dirichlet problem for a plane domain whose 
boundary curve has a tangent and finite curvature at each point. Chapter III deals with symmetric 
kernels K(z, y) = K(y, x), and orthogonal systems of functions. Among other topics, one finds the 
Gram-Schmidt orthogonalization process; the Riesz-Fischer theorem; the bilinear formula for a kernel 
K(z, y) in terms of its eigenfunctions and eigenvalues; the Hilbert-Schmidt expansion theorem and its 
applications; and Mercer’s theorem on positive kernels. In the section on the convergence of the Neu- 
mann series, the reviewer missed A. Weinstein’s simple but little known criterion for convergence, 
(see Sitzungs. Berlin Math. Gess. 26, 1927, p. 168) namely that 
~} 


1 
|A | S | max [ | K(x, y) | dy] , 
z “0 


which compares favorably with C. Neumann’s: 


| >| S [max | K(z, y) |]; 
(z,y) 
and is independent of E. Schmidt’s 
1/2 


1 al 
lA|s l | K*(x, y) dx dy . 


70 70 
Chapter IV is entitled: ‘‘Some types of singular or non-linear integral equations”’, and deals in part with 
theories which are still in the process of formation, as the author states. The author returns here to 
“Cauchy principal value integrals’ and presents many of his own contributions. Of particular interest 
to applied mathematicians is the section devoted to the finite Hilbert transformation and the airfoil 


equation: 


f(z) = if roy 


the integration being over —1 < y < 1. 

In closing, an example of the author’s keen concern for clarity might be mentioned. On page 130, 
after the definition of the jump discontinuity of the Green’s function, G, there is a dagger, referring to 
a footnote at the bottom of the page, which reads as follows: ‘‘Be careful. Some authors (including 
myself in [48]) consider this jump with opposite sign and this implies a change of G into —G.”’ 


J.B. Diaz 


(Continued on p. 106) 
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—NOTES— 


SOLUTION OF NON-LINEAR EQUATIONS 


By ANDRE N. GLEYZAL (U. S. Naval Ordnance Laboratory) 


A standard procedure for solving a set of n equations in n variables: 
f:(x;) = 0, j= By *** pity (1) 


where n is an integer and the f; are real functions of the variables x; may be described 
as follows. We start with a trial point P as a first approximation, and expand the function 
f, in the vicinity of P. Retaining only the linear portions of the expansions we solve 
the resulting linear equations obtaining a point Q. The point P is now replaced by Q 
and the calculation is repeated. It has been shown under mildly restrictive conditions 
that this process yields a rapidly convergent series if the trial point P is sufficiently 
close to the solution. 

We consider the procedure where instead of Q we select a point R on the line PQ 
which gives a minimum value to the sum of the squares of the function f; . 


Suppose P = (£;) is an approximate solution. Then, assuming the functions (1) 
are analytic, let 
; of I a’f F 

» = o(dz;) = i(E;) —_dzr,+= ——— dr, dx eee 

. »D [ fi) + p» oz; ° 2! » én, 0z,; °° + , 
where 

dz; = x; — &; , t,j,k,l=1,-++ ,n. 

Let 


¥= var) = > Ee + >a ax, |. 


, 
Thus y differs from ¢ by terms of order two or higher in dz; . Let 
dx; = du; 


be a solution* of the linear equations: 


t Of — 
filé) + > >. dx; = 0. 


Thus Q = (£; + du;). Now consider the two functions of X: 
¢g = g(A du;) = g(A), 
¥ = Adu;) = YQ). 

For \ = 0: 


o= wg, dg/dyX = dy/dx. 


Received November 4, 1957; revised manuscript received January 16, 1958. 
*See Householder, Principles of numerical analysis, McGraw-Hill, New York, 1953, pp. 81-85. 
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Also, (A) is a parabola which takes on its minimum value, zero, at \ = 1. The function 
g(A) is always non-negative. Hence ¢(A) must reach a relative minimum for a value 
\ = A, , where 0 < XA,, S ~. The quantity \,, may readily be calculated and the cor- 
responding point: 


R = &; + Am du;) 


then used as a starting point for the calculation in place of é; . ‘Thus ¢ must decrease with 
every iteration. Moreover, the iterative method described above will converge in one 
step if Eqs. (1) are linear, and it has been found in many calculations that the method 
converges rapidly if the initial approximation is sufficiently close to the solution. 

If the derivatives 0f;/dx; are not readily found one may approximate them by 
difference quotients Af;/Az,; . In extremely non-linear problems where independent and 
dependent variables become very small or very large in the process of calculation the 
choice of suitable increments Az; may not be obvious. In this case one may take the 
increments Az; to be differences in successive approximations of z; . It has been found 
in many calculations that for sufficiently smooth functions approximately maximum 
numerical accuracy is then obtained for the solution. 

The method described here was successful in solving certain extremely non-linear 
chemical equilibrium equations where independent and dependent variables were un- 
predictably small or large. In other methods—such as the “steepest descent method’’— 
which were tried the rate of convergence decreased rapidly and the method became 
impractical as the solution was approached, whereas in this method the rate of con- 
vergence increased rapidly. Combinations of the steepest descent method with this 
method suggest themselves. For example, let 


z; =§& +rAdu; + ud, , 
where u is a parameter and dv; = —grad ¢. The parameters A and u are then adjusted 


so that ¢ is minimum. 


ON CONVERGENT PERTURBATION EXPANSIONS* 
By RICHARD BELLMAN (The Rand Corporation) 
AnD TOMLINSON FORT (University of South Carolina) 


1. Introduction. In this paper, we wish to consider the Sturm-Liouville equation 
uw + Xf) + eg(a))u = 0, ti) 
u(0) = u(1) = 0, 


and the problem of obtaining power series expansions for the first characteristic value 


and first characteristic function. 
Let A, and u(x) be respectively the first characteristic value and associated charac- 


teristic function of the equation 
u’’ + Af(x)u = 0, 
f(x) (2) 
u(0) = u(1) = O. 





*Received February 24, 1958. 
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Let uo(x) be normalized by the condition that u’(0) = 1. Then it is well known how to 
determine perturbation expansions of the form 
A = No + €A, +h, + oe 
(3) 
uU=Uteaut+eu,t+-::-, 
for the smallest characteristic value of (1) and the characteristic function reducing to 


uo When « — 0. 
In [1], it was shown that in a number of similar situations involving linear operator 


equations of the form 
L(u) + «eM(u) = 0, (4) 
under suitable assumptions it was possible to find a function h(e) with the property 
that the perturbation series 
uU = Up + uh.) + uhle)? + --- (5) 


converged for all ¢« > 0, a case of physical importance. The choice of A(e) depends upon 
the a priort knowledge of the location of certain singularities of u regarded as a function 


of e. 
In [3], in response to a question posed in [2], the following result was established. 


Theorem 1. With reference to the equation, assume that 
(a) f(x) and g(x) arecontinuousfor 0<2z< 1, (6) 
(b) f(z), gz) >a >0, for O<2z <1. 

Then, A, the smallest characteristic root of (1) is an analytic function of € for 


Re («) > b = —Min g(z)/f(z). (7) 


Furthermore, the only singularity of A on the line b + ir, —7 <1 < ~,tsatr = 0. 
Using this result, we shall show how to determine A(e), as above in (5), so as to obtain 
an expansion for A which is convergent for all « > 0. 


2. The choice of h(e). Consider the transformation 
6 = «e/(k + 6), (1) 

where k is a real quantity to be specified in a moment. 
The unit circle, | 6| = 1, is carried into the line 2Re(e) + & = 0in the complex e¢-plane. 
Hence, if 2Re(e) + k > 0, we have! 6| < 1. Let us then choose k to be a positive quantity 


in the interval 


0<k< -—b, (2) 

where b is defined by (1.7). 
Then A(e) = A[ké/(1 — 6)] is analytic for | 6! < 1. Hence, we have an expansion 
A(e) = A[ks/(1 — 8] = >> A, 0", (3) 


n=0 
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valid for | 6 | < 1, or, finally, 
Af) = > Awle/(k + OC], (4) 


valid for « > 0, and, actually, for Re(e) > —k. 


3. The expansion for the associated characteristic function. Consider the solution of 


u’’ + A(e)[f(x) + eg(x)Ju = 0, (1) 
u(0) = 0, u’(0) = 1. 


The classical iteration procedure for obtaining u shows that w is an analytic function of 
A(e) and efor 0 < x < 1, for all € and A(e). Consequently, uv is an analytic function of 


e for Re(e) > 6, uniformly in 0 < x < 1, and thus possesses an expansion 


“uo » u,(x)[e/(k + ©)]”, (2) 
valid for e > 0. 
4. Perturbation procedure. To determine the sequences {\,} and {u,(x)}, we 
employ a standard perturbation technique starting with the equation 
u’’ + Alf(xz) + kég(x)/ — 5)Ju = O, 1) 
( 


u(O) = u(1) = O. 


REFERENCES 
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THE NUMERICAL SOLUTION OF AN INFINITE SET OF LINEAR 
SIMULTANEOUS EQUATIONS* 


By B. NOBLE (The Royal College of Science and Technology, Glasgow, Scotland) 


The usual method for numerical solution of 


>> 4,.2, = b,, (7 = 1,2, -°:), (1) 
& 1 


is to solve the first nm equations in nm unknowns for some fixed n, and to regard the results 
as approximations to the x, . The accuracy of this approximation is invariably doubtful. 
A slight modification of the Choleski’ or Crout® method for solving a set of simultaneous 
linear equations can be used to provide a systematic procedure for solving the first 


m X m equations with m = 1, 2, 3, --- in succession. It is then possible to estimate in 


various ways the error involved in solving only a finite number of equations. The following 
discussion should be intelligible from either the Crout or the Choleski point of view. 
*Received February 14, 1956; revised manuscript received March 7, 1958. 
1L, Fox, H. D. Huskey, and J. H. Wilkinson, Quart. J. Mech. and Appl. Math. 1, 149-173 (1948). 
2P, D. Crout, Trans. Am. Inst. Elec. Engrs. 60, 1235-40 (1941). 
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Using matrix notation, the equations AX = B are solved by decomposing the matrix 

A into a product of triangular matrices A = LU where L is a lower triangular matrix 

and U is an upper triangular matrix. The following systems are then solved in succession: 
LY = B: UX = Y. 

We shall adopt the layout used by Crout in which we start with an “original matrix” 


D and write down an “auxiliary matrix”’ E: 





| A » Qin, °°? » Am, O, |b » U2 » Uin 1% | 
D = | do , Gen, *** 5 Gon y Oy |: E = | ly, ’ be, °°° » Uan » Yo |, 
Ant ’ Ano ne tale ] Onn ; b,, | Ll ls ’ hie ’ a 


where the l,, , u., , y, , are elements of L, U, Y, respectively. For instance U consists 
of zeros below the principal diagonal, units in the diagonal and the u,, above the diagonal. 
Elements of the auxiliary matrix are determined in the following order: 
(i) Elements in the rth column on and below the principal diagonal (i.e. 


rr» : ¥ee 2 


l,r 


(ii) Elements in the rth row to the right of the principal diagonal (i.e. 


u , tee ae So 


This sequence is performed for r = 1, then for r = 2, and so on, up to r = n, by 


means of the formulae 


r—1 


a—1 , 1 
lL, =a,,- © a ee Ure = 7 | Ore — > re ae oe | 
m=1 rr m=1 


The proposed method replaces the solution of UX = Y by the following procedure 


(x,"’ is the mth approximation to z, , etc.). 
(a) Take as a first approximation to the solution of the infinite system: 
(1) (1) 
t, = Y, : z, =0, (ry = 2,3, --:-). 
(b) Obtain a second approximation by solving: 
(2 (2) 
z + U2, = 0 
(2) 
ts: = We - 
Then 
(2) (1 (2) (2) (2) ‘ 
ag” 6 le: ts = ¥,: zs, «= @, (yr = 3,4, ---). 
third approximation by solving: 
(3) (3) (3) 
a + Ut. +323 = 0 


(3) (3) 
Zo + Uses = 0 


(ec) Obtain : 


Q 
- 3) 


fs @* 9s; 
Then 
(3 (2) (3) (3) (2) (3) 
4 = 2; ade le =@=Z. 2: 


(3) (3) _ 
a =y,: 2 = 0, (r = 4,5,6,---) andsoon. 
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The theory behind the method is obvious if we reflect that at the mth stage we have 
obtained the exact solution of the first m equations in m unknowns. 
As an example consider the following original matrix which is derived from the first 
four equations of an infinite system* (the elements in the check column are just the 


sums of the elements in the corresponding rows of the original matrix): 


The auxiliary 


2; Ls Ls Le Constant Check 
+5.665,118 —0.240,000 +0.059,172 —0.022,400 +2.670,644 8.132,534 
—0.240,000 +1.270,836 —0.103,806 +0.049,941 —0.049,299 0.927,672 
+0.059,172 —0.103,806 +0.708,321 —0.051,132 +0.006,400 0.618,955 
—0.022,400 +0.049,941 —0.051,132 +0.489,615 —0.001,666 0.464,358. 


matrix is: 


+5.665,118 —0.042,365 +0.010,445 —0.003,954 +0.471,419 1.435,545 
—0.240,000 +1.260,668 —0.080,353 +0.038,862 +0.050,641 1.009,150 
+0.059,172 —0.101,299 +0.699,563 —0.067,129 —0.023,393 0.909,478 
—0.022,400 +0.048,992 —0.046,961 +0.484,470 +0.010,969 1.010,969. 





The procedure described above then gives 


x, = +0.471,419. 

z, = +0.002,145: 23? = +0.050,641. 

2; = +0.000,165: 2,” = —0.001,880: 2z;” = —0.023,393. 

2 = +0.000,020: z{ = —0.000,367: 2S = +0.000,736: 2z{” = +0.010,969. 


; ; = aa ; 
» agree with those given in* to a unit in the sixth decimal place. 
accuracy can be obtained by comparing the magnitudes of 


The values for z; 
In addition, an estimate of 


the corrections for successive values of m = 1, 2, 3, --- . By inspection of the above 
results it appears that the approximate x, obtained by solving the 4 X 4 set of equations 
are accurate to roughly 5, 4, 3, 3 decimal places for s = 1, 2, 3, 4 respectively. 


When the corrections z,;”’ do not converge rapidly to zero it will be desirable to give 
definite upper and lower limits for the unknowns. This may be easy if the corrections 
alternate in sign. In more difficult cases the following procedure has proved useful. 
Suppose that we have found the estimates x}”’ of x, by solving the first m equations in 
m unknowns. Suppose that from these or otherwise it is possible to estimate the x, for 
s=m+1,m-+ 2, ---. Denote these estimates by X}”’(s > m). Then a second estimate 
of the z, is given by the quantities (2{” satisfy 


> «ai =— >» Oaage 


=] s=m+l1 


a) 


’ + u;”’) where the u;” 
(r = 1 to m). (2) 


For sufficiently large m the u{” are small corrections and a comparatively crude estimate 
of the sum on the right hand side of this equation is adequate. If convenient we can 





8A. Weinstein and D. H. Rock, Quart. Appl. Math. 2, 262-266 (1944). 
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replace the a,, by an approximation with the correct behavior for large s. It will often 
be possible to replace the summation by integration. 
As an example consider the following equations’: 


2.7592, + 0.3552, + 0.3022; + 0.2672, = 1.0000 
0.21427, + 4.457x, + 0.2272, + 0.21427, = 0.3333 
0.1432, + 0.1782, + 5.691z; + 0.1782, = 0.2000 
0.107z, + 0.1422, + 0.1517; + 6.70327, = 0.1429. 


These are the first four equations in four unknowns from an infinite system for which 


b. = (2r — 1)" and it is known that for large r, 
a,, && (2n)'7(2r — 1)'7: ~—a,, FY (28 — 1)'/{(2n)'7r +8 — 1}. (3) 

Also all unknowns are positive so that by inspection 

a, < a,b, = (2n)7'7(2s — 1)7*”. (4) 
We shall use the following estimate in (2): 

Xi” w 2 (2m — 1)°7(28 — 1)7*”. (5) 
This will give an overestimate of the magnitude of the right hand side of (2) since from 
(4) the xz, decrease at least as rapidly as (2s — 1)~*” and from the numerical results 
given below it will appear that the x{” are overestimates of the z,, . We replace the 


right hand side of (2) by 


tw) (2m — 1)” [ dx 
ative (22 — IG +2 — 1) 


™ (Qn) } 
ae (2m — 1)°*” | (2 + 2m — ‘) 
7m (Qm)'7(@r — 1) 8 on ~C 


The final results are presented in the following table where the upper figures give 
xi” | the lower give 2” + u{”. The convergence of the results from above and below 
as m is increased indicates that upper and lower limits have been established for the 
unknowns. But the convergence is much slower than in the example given above and 
it would be difficult to establish accurate estimates of the unknowns from the upper 








limits alone. 


a) 0.3624 
gg = 


1 = 


0.3412 
(2) 0.3550 (2) 0.0577 
xi? = ry” = 
0.3459 0.0531 
(3) 0.3524 (3) 0.0566 (3) 0.0246 
zy = To = T3 = 
0.3469 0.0536 0.0225 
4 0.3512 (4) 0.0560 (4) 0.0242 (4) 0.0140 
“4 = —~= = = = “4 = 
0.3474 0.0539 0.0226 0.0128. 





4D. S. Jones, Proc. Roy. Soc. 217, 153-175 (1953). 
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We next consider evaluation of the following infinite series (see*): 
§ = > kz, , 
o=1 
where 
ky* = 2.33, = 4.22, kk’ =65.51, kz’ = 6.55, 


1/2/< 1/2 ¢ 
(2s — 1)” for large s. If we evaluate 


S, = > k.a\” 


s=1 


and k>* ~ (27) 


we obtain the following results: 
S, = 0.155, : S, = 0.166, : S; = 0.169, : S, = 0.170; , 
These underestimate S. In order to obtain overestimates of S we evaluate, using the 
previous approximation (5), 
. ai 2 < 1 
Si = > ke” + k,x(2n — 1) z. (2s — 1? 
s=1 _ 


s=n+l1 \ 
, ( n) /< 2/ -1 
~ ; ka.” + k,x,"'(2n — 1)*(4n)™, 
s=1 
where we have replaced summation by integration as before. This gives 
S; = 0.194, : S3 = 0.181; : S3 = 0.178; : Si = ©.176, . 


The mean values 3(S, + S/) are remarkably constant. We can say from these results 
that 0.171 < S < 0.177 with a probable value of S = 0.174. 


ON SIMPLE SUBHARMONICS* 
By C. 8. HSU** (University of Toledo, Toledo, Ohio) 


Introduction. In a recent paper [1], Rosenberg clarified considerably the curious 
subharmonic phenomena of non-linear oscillations by introducing the concepts of strong 
subharmonic solutions and of the simple subharmonics. He considered a system with a 
single degree of freedom, whose mechanical model might be a mass under the action of 
an elastic force, linear or non-linear, and of a simple harmonic forcing function of fre- 
quencey w. The equation of motion of the system can be put in the form 

d’x 


Te + f(x) = P» cos at, (1) 


where P, is proportional to the amplitude of the external forcing function. 





*Received June 16, 1958. 
**Now affiliated with University of California, Berkeley. 
+Numbers in square brackets refer to the bibliography at the end of the paper. 
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If we are concerned only with the steady periodic solutions of (1), we can express 
them in terms of Fourier expansions. If such an expansion contains a term A, cos 
(w,t + ¢,), where w, = w/r and r is an integer, and if A, + 0, the solution is said to be 
subharmonic of order 1/r, and we have 


x = A, cos (2 t+ e,) + > A; cos wt + ¢)). (2) 
tr 
Rosenberg called the solution a strong subharmonic of order 1/r, if | A, | >> | A, | for all i. 
He called the solution a pure subharmonic if A, ~ 0 and A, = 0 for all z. A solution is 
said to be a simple subharmonic if it is a pure subharmonic with zero phase angle ¢, , 
ie. Z = 2% cos (wt/r). 
tosenberg showed that, for a class of differential equations of the following type 


r 


a +rvt+k >> alt" = k cost (ry = 1, 2,3, ---), (3) 


periodic solutions are of the simple subharmonic type 
£ = cos 7/r. (4) 


Here a set of new variables has been introduced. In terms of the original variables 
they are: r = wt, £ = x/z) and k = P,/(x w’). A table was included in [1] giving the 
values of a,” in (3). The linear problem is a special case of (3) and (4). Then Eqs. (3) 
and (4) respectively become 


© + (1 + WE = k 008 7, 6) 
dr 


and 


— = cos fr. (6) 


This is the usual linear response solution. 

In the present note the writer will try to amplify and discuss some aspects of this 
interesting subharmonic problem. It will be shown that by a single transformation all 
of the equations (3) and their subharmonic solutions (4) can be derived from the linear 
case (5) and (6). 

Tchebycheff polynomials. A careful study of Table 1 of [1] which gave the values 
of a” indicates that a{” are simply the coefficients of various Tchebycheff polynomials 
of the first kind. This was also mentioned in another paper [2] by Rosenberg. Thus 
Eq. (8) can be written as 


£54 t+ ee «bees (7) 
dr 


and its periodic solution given by the simple subharmonic (4). For ease of reference 
the Tchebycheff polynomials of the first kind will be briefly described here. The simplest 
definition is as follows. If we put y = cos 6, the Tchebycheff polynomial of the first 
kind of order r, denoted by T,(y), is defined as 


T.(y) = cosré. (8) 
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The first few of these are: 


Ty) = 1, T;(y) = 4y’ — 3y, 
Ti(y) = y, T.(y) = 8y* — 8y + 1, (9) 
T'.(y) = 2y” — l, CCT eee ee ee Tee 


We note that if § = cos (r/r), the sum of the first two terms of (7) is identically 
zero and the third term is precisely k cos 7, from the definition of 7’, . Thus (7) is satisfied. 
Physically, this means that with a properly constituted elastic non-linearity there exists 
an amplitude of the forcing term such that a simple subharmonic response will cause the 
spring force to balance partly the inertia force and partly the external force all the time. 

It is easily seen that (7) is not the only differential equation which possesses simple 
subharmonic periodic solutions. As a matter of fact we have the equation 


dt 


qe +r + k(—1)'"T,(¢) = k cost, (r even integers) (10) 
a 


with periodic simple subharmonic solutions 


t = sin—, (11) 
r 


and the equation 


- +re + k-1)°°"°T (6) = ksin + (r odd integers) (12) 
* 


with solutions also given by (11). 

Transformation with Tchebycheff polynomials. Let us now take a different approach 
to arrive at Eq (7). Let us return to the linear problem with a forcing term, i.e. Eqs. 
(5) and (6). If we apply the following transformation to (5) and (6) 


§=T,(n), (13) 
where 7 is a function of r, Eq. (5) becomes 


dr- 


+ T.(n) + kT.(n) = k cost (14) 


and the solution is, from (6) and (13), 


n = cos ns (15) 


Equation (14) may be rewritten as 
oo d’n 1 d’n lm - | - ; 
| T | bes E + r n| + dr + r -T kT (n) = k cos T) (16) 


by adding and subtracting the terms in the second square bracket. The first two square 
brackets are, however, identically zero with ¢ and n given by (6) and (15). It follows 
that » = cos (7/r) must be a solution of 
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d’n 


dr: 


+r-?n + kT An) = k cos 7, (17) 


which is the same as (7). 


Similarly, wutting € = (— 1)” 7,(n) in (5) we can arrive at (10) and (11). Equa- 
tion (12) ana ..s solution can be established in the same manner by using a slightly 
different form of (5). 

This approach shows that by a single transformation all the non-linear differential 
equations possessing simple subharmonics as their solutions can be derived from the 
linear problem which we understand much better. In other words, the simple sub- 
harmonic solutions might be considered merely as different representations of the same 
linear response solution in various frameworks. The physical significance, if any, of 
this transformation is not apparent at the present time. 

Uniqueness of the solutions. There has been some discussion among interested people 
concerning the uniqueness of the simple subharmonics as steady-state solutions of (7). 
Since it is known that the linear response (6) is a unique steady state solution of (5), 
and since Eq. (7) and its simple subharmonic solutions can be derived from the linear 
case by a proper transformation, which is one-to-one when continuity conditions are 
imposed on the solutions, it may be concluded that these simple subharmonic solutions 
are unique so far as the steady state is concerned. 
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Correction to my paper 


ON THE DAMPED OSCILLATIONS EQUATION WITH 
VARIABLE COEFFICIENTS 


Quarterly of Applied Mathematics, X VI, 90-93 (1958) 
By E. V. LAITONE (University of California, Berkeley) 
Equation (13) should have been written 
$(0) > {o(O)u(O)? + fu’(O) + p(O)u(O)/2]*} (13) 


that is, the u’(0) should not be squared. Similarly the last equation on page 93 should 
contain only a’(0) in the brackets. 
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BOOK REVIEWS 
(Continued from p. 94) 


An introduction to the theory of random signals and noise. By Wilbur B. Davenport, Jr. 
and William L. Root. MeGraw-Hill Book Co., Ine., New York, Toronto, London, 
1958. ix + 393 pp. $10.00. 


Let us begin by complimenting the authors upon having written a very clear and readable intro- 
duction to the mathematical theory of random signals and noise. By writing in a leisurely and lucid 
fashion, and by means of a host of well chosen examples, they provide an almost painless entry into one 
of the most significant fields of modern analysis. 

A survey of the chapters will indicate the scope of the book. 

The first chapter is introductory. The next five chapters provide the background of probability 
theory required by the reader, elementary aspects, random variables, probability distributions, random 
processes and statistical averages. There is a discussion of sampling and the central limit theorem, 
followed by a chapter devoted to spectral analysis. 

The next chapter is devoted to a discussion of shot noise in vacuum tubes and contains a detailed 
analysis of the physical origin of many of the problems treated in the book. This is a valuable addition 
to the mathematical treatment. 

The following two chapters are devoted to Gaussian processes and linear systems. The problem of 
determining optimal linear filters leads to a discussion of the Wiener-Kolmogoroff theory of statistical 
detection of signals. 

The only objection of any consequence that we wish to make is to the use of the “impulse function.” 
If the same notation is desired, it is best perhaps to refer to the distribution functions of L. Schwartz. 


Otherwise, the Riemann-Stieltjes integral can be used. 
RicHARD BELLMAN 


Theorie der Beugung Elektromagnetischer Wellen. By W. Vranz. Springer-Verlag, Berlin- 
Géttingen-Heidelberg, 1957. 123 pp. $5.20. 


This monograph, No. 4 of the series ‘Ergebnisse der angewandten Mathematik’’, gives a compre- 
hensive treatment of the solution of Maxwell’s equations for the boundary conditions appropriate to 
diffraction problems. It is divided into three parts. The first part discusses general mathematical ap- 
proaches, introduces the Green’s dyadic as analog of the Green’s function of the scalar field, and formu- 
lates the needed boundary conditions, The second part treats diffraction by bodies with smooth surfaces, 
deriving the exact solutions for the cylinder and sphere in detail and indicating the application of the 
integral-equation method to diffraction at a weakly curved surface of arbitrary shape. The last part 
deals with diffracting bodies having sharp edges. Sommerfeld’s exact solutions for the perfectly conduct- 
ing wedge and halfplane are derived and approximate methods for the treatment of diffraction at slits, 
apertures, and circular disks are discussed. A full bibliography aids the reader in locating approaches 
to the diffraction problem not covered by the text. The book constitutes a valuable addition to the 
literature on the mathematical theory of diffraction by an author who has contributed liberally to the 


subject. 
Ik}. G. RAMBERG 


Mechanical resolution of linguistic problems. By Andrew D. Booth, L. Bandwood, and 
J. P. Cleave. Academic Press, Inc., New York, and Butterworths Scientific Publica- 
tions, London, 1958. vil + 306 pp. $9.80. 


This volume, the senior author of which is one of the pioneers in the field of machine translation, 
is, largely, an account of experiments at the Birkbeck College Computation Laboratory in the application 
of digital calculators to linguistic problems. The machine used was the APEXC all-purpose electronic 


computer which is a laboratory model of the somewhat more versatile HEC General Purpose Electronic 
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Computer, a business machine produced by the British Tabulating Machine Company. The input on 
the APE-XC is in the form of a punched paper tape, and the linguistic information is encoded in binary 
form, with the international teletype code utilized for this purpose. The authors recognized that the 
teletype code is not the most suitable one for application to linguistic analysis but adhered to it for 
practical reasons, mainly because ordinary teletype equipment could be attached to the machine in 
such a case, and the input and output procedures thus facilitated. The storage capacity of the computer, 
which takes the form of a magnetic drum, is 8,192 words of 32 binary digits each. The limitation of 
word length to 32 binary digits, i.e. to no more than six letter words represented in the five digit binary 
code, means, in practice, that more than one storage location must be utilized for the many longer 
words of a language. The APEXC ean perform basically the same operations as the other standard 
computers. 

The book contains twelve chapters and subject and name indices. It lacks, unfortunately, an anno- 
tated bibliography, although footnotes to chapter one, which is a brief history of computer application 
to linguistic problems, contain references to the most important works on the subject. Chapter two gives 
some basic information about the nature of calculating and data processing machines and the following 
chapters turn to the specific linguistic problems in which such machines may find suitable application. 
The authors point out that there are, besides translations, many other important problems of linguistic 
analysis for which the computer may be utilized to advantage. Some of these problems and their pro- 
gramming are discussed: the word counts and concordances, the analysis of sentence structure, problems 
of lexicography and syntactical investigation, and literary dating such as, for example, the Platonic 
chronology. In chapter six, the authors turn their attention to mechanical translation and its program- 
ming. The complications which may arise in such analyses, such as the treatment of stems and endings, 
compound words, multiple meanings, or idioms, are discussed but not always thoroughly explored. 
Chapter seven examines the possibility of mechanical transcription of Grade II Braille (the abbreviated 
form which uses special symbols for frequently occurring groupings of letters). Technically, such trans- 
cription is found feasible in spite of some ambiguities, but, because of the low output speed, it appears 
impractical economically. Chapters eight and nine suggest experimental schemes for the translation of 
French and German respectively, while the following chapter gives a brief account of the work of the 
Academy of Sciences of the USSR on mechanical translation of English into Russian. The question of 
multilingual translations is very briefly mentioned in chapter 11; the final chapter contains suggestions 
of technical details of a machine which would be designed specifically for translations. It emphasizes 
that some operations of presently available computers would be unnecessary (multiplications), while 
others, now lacking, should be added, such as greater word length capacity. The all-important question 
of adequate storage capacity is duly stressed, and a table of a proposed instruction code for such a 
linguistic machine is included. 

The main value of this book, in the opinion of this reviewer, lies in the clarity with which it summarizes 
some of the basic problems concerning the application of digital calculators to linguistic problems and 
shows, in proper perspective, the limitations of the presently available machines. However, not all 
parts of the book will be equally informative for all potential readers. The linguist who may be drawn 
to the book will often be left puzzled by the lack of clarity of the programming descriptions, and perhaps 
somewhat skeptical of the way in which difficult linguistic problems have been glossed over. The analysis 
of French and German would have profited from the application of some pertinent principles and tech- 
niques of structural linguistics. Neither will the reader interested in details of programming on the 
APEXC find it in the present volume; instead, he is referred to another book which offers this information. 


Henry Kucera 


Proceedings of the Third Midwestern Conference on Solid Mechanics. University of Michigan 

Press, 1957. vi + 250 pp. $5.50. 

The book under review contains the following 15 papers read at the Third Midwestern Conference 
on Solid Mechanics, held at the University of Michigan in April, 1957. 

1. Pitching Instability of Rigid Lifting Surfaces on Viscoelastic Supports in Subsonic or Supersonic 
Potential Flow, by H. H. Hilton; the investigation of single-degree-of-freedom torsional instability of 
rigid lifting surfaces on flexible supports. 2. On the Motion of Thermo-Viscoelastic Solid, by M. Lessen; 
equations describing small motions of a viscoelastic solid are set down from a thermodynamic point of 
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view and the propagation of longitudinal waves is studied. 3. Approximate Yield Conditions in Dynamic 
Plasticity, by P. G. Hodge and B. Paul; a simply supported circular cylindrical shell is subjected to a 
blast-type radial pressure loading, and the resulting plastic displacement is computed according to 
square and hexagonal interaction curves. 4. Damping of Vibrations in Elastic Rods and Sandwich Struc- 
tures by Incorporation of Additional Viscoelastic Material, by H. I. Plass; the paper is concerned with a 
method of reducing the vibrational energy transmitted through a structural member (circular bar with 
viscoelastic layer, and sandwich panel having viscoelastic core). 5. Thermal Stresses in an Infinite Plate 
of Arbitrary Thickness, by E. L. McDowell; the problem of determination of the steady-state thermal 
stresses and displacements in an isotropic elastic region bounded by two parallel infinite planes is studied. 
6. Elastic Stability of Conical Shells Loaded by Uniform External Pressure, by C. E. Taylor; an approxi- 
mate solution of large deflection theory for buckling pressure for complete conical shells is obtained by 
Rayleigh-Ritz method. 7. Torsion with Warping Restraint of Tapered Beams, by M. Goulard, Hsu Lo 
and R. I. H. Bollard; an exact and an approximate solution for a cantilever thin-walled tube of rec- 
tangular cross section in the shape of truncated cone is derived. 8. A Theoretical and Experimental Investi- 
gation of the Stability of a Hovering Helicopter Rotor Blades, by J. Zvara; model tests and a consistent 
theoretical flutter analysis of a system having three degrees of freedom are made. 9. Reflection and 
Transmission of Elastic Pulses in a Bar at a Discontinuity in Cross Section, by E. A. Ripperger and H. 
Norman Abramson; reflection and transmission coefficients based upon an eler-entary theory and steady- 
state conditions are compared with experimental results. 10. The Timoshenko ?-am on an Elastic Founda- 
tion, by S. H. Crandall; a dynamic response of Timoshenko beam supported by a massless elastic founda- 
tion to a fixed pulsating load and to a moving constant load is obtained. 11. On Rotating Blades with 
Flexible Mounting, by Hsu Lo and C. E. Danforth; static and free vibration problems of such systems 
are investigated. 12. The Bending Vibrations of a Twisted Rotating Beam, by W. P. Targoff; an analytic 
method is developed which permits the computation of the natural modes and frequencies of rotating 
beams, the results of computations are compared with experimental results. 13. Effect of Damping on 
Vibration Frequencies of Simple Systems, by L. E.. Malvern; one-mass force vibrating system with viscous 
damping across part of spring or torsional pendulum with friction bearing are studied. 14. Bending- 
Torsion Flutter Sensitivity in Incompressible and Supersonic Flow, by A. S. Richardson; it is shown that 
“flutter sensitivity’’ (the slope of the velocity vs. structural damping curve at the flutter point) is not 
an exact measure of the violence (or lack of violenee) of the flutter mode. 15. Flutter of Curved Plates 
with Edge Compression in a Supersonic Flow, by Y. C. Fung; the importance of initial warping and pres- 
sure differential across the plate on the panel flutter characteristics is demonstrated. 

The papers contained in this book, as a rule, are of significant theoretical and applied importance. 
It is to be regretted that the authors are not always aware of the corresponding results published in 
Russian. For instance, E. L. McDowell seems to be unaware of the results obtained by G. N. Maslov 
(An Elastic Problem of Thermoelastic Equilibrium. Izvestia Vsesouznogo Instituta Hydrotechniki, 
v.23 (1938)) and C. E. Taylor seems to be unaware of the results of E. I. Grigoljuk (Nonlinear Vibrations 
and Stability of Shallow Columns and Shells. “Izvestia’”’ of Academy of Sciences, USSR, Department 
of Technical Sciences, No. 3, 33-68 (1955); Large Deflection Instability of a Closed Sandwich Conic 
Shell under Uniform Normal Pressure. Inzhenernyi shornik, v.22, 111-119 (1955)). 


G. 5. SHAPIRO 


The Scientific Papers of Sir Geoffrey Ingram Taylor. Volume I: Mechanics of Solids. 
Edited by G. K. Batchelor. The University Press, Cambridge, 1958. x + 593 pp. 
$14.50. 


The present volume contains forty-one papers on a wide range of topics in mechanics of solids. 
They are ordered chronologically starting with the well-known paper on the use of soap films in solving 
torsion problems, which was written in collaboration with A. A. Griffith and published in the Proceedings 
of the Institution of Mechanical Engineers in 1917, and ending with a paper on strains in crystalline 
aggregates that was presented at the Colloquium on Deformation and Flow of Solids in Madrid in 1955. 
Some papers written for Government agencies are published for the first time in this volume. These 
include: Notes on the ‘‘Navier Effect’’ (1925); Lattice Distortion and Latent Heat of Cold Work in 
Copper (1935); Propagation of Earth Waves from an Explosion (1940); Calculation of Stress Distribution 
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in an Autofrettaged Tube from Measurement of Stress Rings (1941); The Plastic Wave in a Wire Ex- 
tended by an Impact Load (1942); The Mechanical Properties of Cordite during Impact Stressing 
(with R. M. Davies, 1942); and The Distortion under Pressure of an Elliptic Diaphragm which is Clamped 


along its Xdge (1942). 


Queues, inventories and maintenance. By Philip M. Morse. John Wiley & Sons, Inc., 
New York, and Chapman & Hall, Ltd., London, 1958. ix + 202 pp. $6.50. 


When a system involves an ordered process moving from stage to stage it is usually spoken of as 
involving flow. In the widest sense flow problems arise in all branches of industry; in the supply of raw 
materials and labour, on the production lines, in the administrative handling of documents such as 
accounts, in the smooth provision of information to management and in the development of the business 
organisation itself. The flow may concern itself with physical quantities such as water behind dams or 
fuel to a factory; or with services such as the issue of air passages, the repair of vehicles, and the pro- 
vision of unoccupied telephone lines. 

One of the standard problems in this field is to ensure that capacity can handle the average flow 
without causing obstruction. To have too much capacity is wasteful; to have too little frustrates the 
process. Somewhere there lies an optimal point which obviates serious risk of bottlenecks without undue 
expense. The problem is to find it. The system may be subject to random disturbance or, more generally, 
to stochastic elements, as when customers arrive at a service counter haphazardly or variations in weather 
may create sudden changes in demand for a product. Thus the flow control often incorporates an element 
of probability theory. The dynamic plus the stochastic element constitute the basic features of the class 
of studies known variously as “queuing theory,’’ “stock control,’’ “inventory problems’’ and so forth. 
They may be considered as a fundamentally important part of operational research. 

Dr. Morse has written an excellent text on the subject. Chapter 1 deals with representation in terms 
of probabilities of state. This leads on naturally to frequency distributions of arrival time at a service 
point and of service times themselves. Chapters 3 and 4 then deal with single and multiple exponential 
channels (the exponential case being one of the few which can be handled with comparatively simple 
mathematics). Chapter 5 then considers the simulation of more complicated cases. There follow four 
chapters on general queuing theory. The text closes with a chapter of application to inventory control 
and one on the maintenance of equipment. There follow a glossary, some tables and a bibliography. 

Anyone who wants an introduction to queuing theory or to operational research can hardly do 
better than begin with this book. The basic ideas are simply and clearly expounded with a sufficient but 
unpretentious mathematical content. It began, so the author tells us, as an introduction to the tables 
with which the work concludes, but rapidly outgrew the prefatory stage and has become a monograph 
in its own right. Altogether it is an admirable piece of exposition in a domain where good expository 
accounts are rather rare. Dr. Morse promises further volumes on computational methods and detailed 
solutions of some of the problems touched on in this book, and if he can maintain his present standard 
he will have made a notable addition to the literature of operational research for which we shall all be 
grateful. 

M. G. KENDALL 


Numerical analysis. By Kaiser 8. Kunz. McGraw-Hill Book Co., Ine., New York, 
Toronto, and London, 1957. xv + 381 pp. $8.00. 


This is an elementary treatise, written principally for engineers, although there is an apparent and 
laudable attention to justification of the methods presented. Nearly half the book is devoted to inter- 
polation, differentiation and integration, including chapters on summation of series and multivariate 
interpolation not usually found in textbooks; much use is made of the lozenge diagram in a manner due 
to the author. There is only cursory mention of the important field of Gaussian quadrature. 

Two chapters deal, respectively, with starting and continuing the solution of ordinary differential 
equations: the methods of Taylor, Picard, Euler and Runge-Kutta serve to illustrate the former, that of 
Milne, principally, the latter. There is attention to the accumulation of errors. Since only fifty pages 
suffice for these two chapters, the selection of topics had to be drastic. 
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A chapter is devoted to the solution of simultaneous equations and the evaluation of determinants; 
there is no attempt to deal with eigenvalue problems either of systems of algebraic or of differential 
equations. The sections on elliptic and parabolic partial differential equations manage to impart a 
surprising amount of information in restricted space, but the chapter on hyperbolic equations is perhaps 
misnamed since the method of characteristics is not mentioned and only a simple problem of the wave- 
equation is solved by finite differences. Fredholm and Volterra equations are discussed in the last chapter 
on integral equations. 

Extensive references to sources are given in footnotes. The work abounds in useful worked and 
unworked examples, is very clearly written and should be easy reading for anyone with a knowledge 
of calculus and differential equations. 

Wa ter F., FrReErBERGER 


Analysis and control of nonlinear systems. By Y. H. Ku. The Ronald Press Co., New 
York, 1958. vii + 360 pp. $10.00. 


This book is a book about engineering problems written by an engineer for engineers. As such, it 
possesses the anticipated virtues and the expected faults. 

On the positive side, it contains very detailed discussions of a number of important examples, 
together with graphs and tables. In addition, it is very clearly written with an extensive bibliography. 

The author constantly emphasizes the necessity for dealing with nonlinear equations and the improve- 
ment in performance that can be obtained by using nonlinear properties. 

On the negative side, there is the objection that the book is too much like a cookbook. Only occa- 
sionally is any effort devoted to an explanation and motivation of the mathematical techniques that are 
employed. Some of these, in addition, are rather clumsy. Thus, for example, the proof of the superposition 
principle given on page 7 makes formidable a property which is an immediate consequence of linearity 
of the equation. The discussion of the Blasius equation on page 206 uses a very crude method of successive 
approximations whereas the method presented by Weyl, Proc. Nail. Acad. Sci., 1941, is much simpler 
and more elegant, and converges rapidly and rigorously. On page 253, the same notation is used for a 


function and its Laplace transform. 


Despite the preceding r¢ marks, we feel that the book could be quite useful if used in conjunction 
with a course given by an experienced instructor. 
RicHarD BuLLMAN 
Studies in the mathematical theory of inventory and production. By K. J. Arrow, S. Karlin, 
and H. Scarf. Stanford University Press, Stanford, California, 1958. x 340 pp. 


$8.75. 


Two of the most interesting classes of dynamic programming processes, viewed from the vantage 
points of both analysis and application, are those of inventory control and production smoothing. 

The typical inventory control process may be described in the following general terms. We possess 
various quantities of different items for which there is a demand of stochastic nature from time to time. 
Since there is a penalty of some type attached to not being able to satisfy this demand, at various stages 
additional quantities of these items are ordered, at costs dependent upon types and quantities ordered, 
the times of ordering, the rate of delivery desired, and other factors as well. The problem is to determine 
the ordering policies which are optimal with respect to preassigned measures of efficiency. 

Production smoothing processes are of quite similar structure, with the additional feature that there 
may or may not be costs attached to varying rates of production and rates of ordering. 

30th are representative of the kinds of multi-stage decision processes that arise in economic and 
industrial activities. 

In order to assess the role of the present volume in the continuing research in these fields, let us 
briefly review the very recent history of the optimal inventory problem. In the first place, there is the 
recognition of this decision process as a problem worthy of attention for its own sake, apart from the 
larger activity within which it is always imbedded. For a history of this phase, see the book by T. Whitin, 
The Theory of Inventory Management, Princeton University Press, 1953, which also contains a discussion 


of various preliminary mathematical attacks upon particular aspects of the general problem. 
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The next phase involves the analytical formulation using the techniques of the modern approach 
to the analysis of stochastic processes. The classic paper of Arrow, Harris and Marschak, “Optimal 
Inventory Policy,’’ Econometrica, vol. XIX, 1951, pp. 250-272, inaugurates the use of functional equations 
to describe the process. 

From this point the road forks. To begin with, there are essential questions pertaining to the existence 
and uniqueness of solution of the nonlinear functional equations obtained in this fashion. These were 
first treated in the paper by Dvoretzky, Kiefer and Wolfowitz, “The Inventory Problem,’’ Econometrica, 
vol. XX, 1952, pp. 187-222, 450-466. Only after these points are settled, can one be sure of a one-to-one 
correspondence between the optimal policies of the original process and the solutions of the defining 
equations. 

Secondly, there is the basic problem of studying not only optimal policies, but also sub-optimal 
policies. Within a class of simple feasible policies, we wish to determine the most efficient. This road is 
followed by Arrow, Harris and Marschak in the paper referred to above, and the inverse problem was 
studied by Dvoretzky, Kiefer and Wolfowitz. 

Thirdly, the functional equation can be utilized to determine the structure of optimal policies as a 
function of the structure of the cost and penalty functions, and also of the probability distribution 
occurring. This study was initiated by Bellman, Glicksberg and Gross in their paper, ‘On the Optimal 
Inventory Equation,’’ Management Science, vol. II, 1955, pp. 83-104; see also R. Bellman, Dynamic 
Programming, Princeton University Press, 1957, Chapter 5. 

In the book under review we find a continuation of this effort, carried out quite elegantly and in 
great detail, and, in addition, a continuation of the description of the outcomes of various approximate 
policies. This last requires the use of the tools of recent renewal and queuing theory. 

One of the most interesting chapters is that devoted to time-lag processes, which is to say those in 
which there is an appreciable delay between ordering and delivery. Uuder the assumption that backlogs 
occur, it is shown that the number of state variables required to describe the process may be drastically 
reduced. Consequently, the results obtained by Bellman, Glicksberg and Gross in the paper cited above 
may be extended to this case. This furnishes a ready approximation to the solution of more complicated 
prot esses, 

Turning to the study of production processes over time, a pioneering work is that of Pierre Massé, 

Les Reserves et la Régulation de l’Avenir dans la Vie Economique,”’ Paris, Hermann, 1946, Volume I, 
Avenir Determimé, Volume II, Avenir Aléatoire. A number of interesting processes give rise via simplified 
mathematical models to variational problems of the following type: ‘‘Maximize a functional of the form 
J(u) = $7, gu, w’, t)dt, subject to combinations of constraints of the type a; < u(t) < a2, a3 < w(t) < ay, 
| 6 h(u, u') dt < a;.”’ The presence of the constraints renders some of these problems of far greater diffi- 
culty than might be supposed upon first observation. Let us note in passing that questions of this nature 
arise with great frequency in the study of current engineering control processes. 

The present series of studies contains an analytic treatment of some producticn processes and of 
some related problems arising in the consideration of hydroelectric systems. References are given to 
other work by Hohn, Modigliani and Simon, by Koopmans, by Karush and Vaszonyi, and by Bellman, 
Glicksberg and Gross. A good deal of unpublished work on these problems was carried out at RAND 
from 1948 on by Paxson, Danskin, and Fleming, stimulated by the book of Massé referred to above, 
Further study of variational problems with constraints may be found in Bellman, Glicksberg and Gross, 
“On some Variational Problems Occurring in the Theory of Dynamic Programming,’’ Rendiconti del 
Circolo Matematico di Palermo, Serie II, Tomo III, 1954, pp. 1-35, and in Bellman, Fleming and Widder, 
‘Variational Problems with Constraints,’’ Annali di Matematica, Serie 1V, Tomo XLI, 1956, pp. 301-323. 

The principal authors, as well as Beckmann, Gessford and Muth, have maintained a high mathe- 
matical level with no sacrifice of clarity. All of the studies are worthwhile starting points for further 
research in this fundamental domain, and the introductory expository section is particularly valuable 
in enabling the reader to understand the origin of the following analysis. 

There are, however, some ways in which the level of the book could be raised even higher. In the 
first place, it would have been worthwhile to have included some computational results indicating the 
dependence of optimal policies, minimum costs and maximum profits upon the different parameters 
entering into the mathematical models. Along these lines, it might also have been helpful to point out 
that the functional equation technique of dynamic programming furnishes computational solutions to 
the variational problems described in the preceding sections in a quite simple fashion without the aid 
of any assumptions concerning convexity, concavity and so forth; see Chapter 9 of the book referred 


to previously, 
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Secondly, it should be emphasized that one of the most important aspects of economic control 
processes of this type is that of predicting the future demand on the basis of the observed past demand. 
We thus encounter “learning”’ or “‘correlative’’ processes of the type that have come into prominence 
with the development of sequential analysis by Wald, and by Wolfowitz, Blackwell, Girshick, Karlin, 
and others. This fundamental problem was attacked in relation to inventory processes in the second of 
the papers by Dvoretzky, Kiefer and Wolfowitz referred to above. 

Considering the size of the book and the amount of material it contains, it is understandable, however, 
that the authors felt that a line must be drawn somewhere. Perhaps, we can look forward to a second 
set of studies devoted exclusively to these deeper problems. In any case, the book is certainly to be 
recommended to anyone working in the fields of mathematical economics and operations research, 


RIcHARD BELLMAN 


Vector analysis. By Louis Brand. John Wiley & Sons, Inc., New York, and Chapman & 
Hall, Ltd., London, 1957. xiii + 282 pp. $6.00. 


This book covers some of the material of the author’s Vector and Tensor Analysis, but at a more 


elementary level. 
L. M. Mitne-THOMSON 





FOURTH U. S. CONGRESS OF APPLIED MECHANICS—1962 


Preliminary Announcement 


The Fourth U. S. National Congress of Applied Mechanics will be held on the 


Berkeley campus of the University of California during June 18-21, 1962. 


Research workers in the theoretical and applied mechanics of solids and fluids 
are invited to submit papers for consideration by the Editorial Committee. Fur- 
ther announcements concerning the preparation of papers and deadlines for sub- 


mission will be made as the Congress draws nearer. 


The members of the organizing committee on the Berkeley campus are: Pro- 
fessor W. GoLpsMiTH, Secretary; Professor E. V. Larronz, Treasurer; Professor 
R. M. Rosenserc, Chairman of the Editorial Committee; Professor W. W. 


Soroka, General Chairman. 
Inquiries regarding the Congress should be addressed to Professor W. Gold- 
smith, Secretary, Division of Mechanics and Design, University of California, 


Berkeley 4, California. 


























